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'  ■  Ray  equations  based  on  an  acoustic  Hamiltonian,  and  formulated  by  Jones  et  al. 
(1986)  for  spherical  coordinates  in  NOAA’s  3-D  ray  tracer  HARPO,  are  adapted  to 
ellipsoidal  coordinates  in  the  oceanic  waveguide.  The  ensuing  modified  HARPO  is 
used  to  model  very  long  range  (up  to  antipodal)  acoustic  paths  in  which  the 
difference  between  geodesics  and  great  circles  is  measurable.  The  eventual  objective 
of  this  modeling  is  to  extract  the  predictable  part  of  the  travel-time  trend  and  fluctua¬ 
tions  along  several  long  paths  that  will  be  used  to  monitor  hypothetical  global  warm¬ 
ing  cffecLs.  The  requirements  for  easy  assimilation  and  representation  of  realistic 
environmental  inputs  are  discussed.  These  requirements,  when  coupled  with  the  pos¬ 
sibility  of  classical  chaos  in  the  ray  paths,  dictate  a  new  software  architecture.  We 
use  the  existing  software,  however,  to  breadboard  and  test  features  of  new  ray  tracers 
in  the  global  boundary  layer,  and  to  support  the  experimental  design  of  a  forthcoming 
pilot  experiment  that  will  use  a  transmitter  located  at  Heard  Island  in  the  Indian 
Ocean  near  Antarctica.  A  path  of  particular  interest  links  Heard  Island,  through  the 
Tasman  gap,  to  the  Washington-Oregon  Coast  in  the  northeast  Pacific.  A  3-D  sound 
speed  model  is  formulated  for  the  Indian/Pacific  Ocean  in  the  region  of  this  path,  and 
a  3-D  antarctic  circumpolar  current  model  with  180  Sverdrup  transport  is  specified  to 
cross  the  path.  Eigenrays  are  computed  to  within  2  m  error  on  18  Mm  in  the  presence 
and  in  the  absence  of  currents.  We  draw  conclusions  specific  to  the  acoustic  com¬ 
munication  channel  on  this  path.  In  particular,  we  infer  a  significant  out-of- 
(gcodc.sic)-planc  deflection  (of  about  5(X)km)  of  the  acoustic  eigenpaths.  This  3-D 
deflection  is  due  to  the  thermohaline  structure  and  can  be  quantified,  at  the  present 
lime,  only  by  ray  tracing. 
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1.  Introduction 

An  acoustic  feasibility  experiment  for 
long-term  monitoring  of  global  warming  trends 
is  planned  for  January  1991.  The  concept  for 
this  remote  sensing  project,  related  to  very  long 
range  acoustic  tomography,  is  described  by 
Munk  and  Forbes  (1989)  and  is  summarized  by 
Gibbons  (1990).  The  single  transmitter  will  be 
at  Heard  Island  (53.5°S,  73.5°E)  in  the 
southwestern  Indian  Ocean.  The  listening 
receivers  will  be  located  in  several  oceans  and 
at  several  distances  (at  megameter  ranges)  from 
the  transmitter.  Windows  for  three  of  the  long¬ 
est  paths  are  westward  from  Heard  Island  to  the 
Atlantic  Coast  of  North  America  and  eastward 
from  Heard  Island  to  the  US  West  Coast.  The 
(unrefracted)  geodesics  that  bound  these  three 
windows  are  shown  in  Figs.  1  and  2.  The  west¬ 
ward,  "Atlantic,"  window  has  the  widest  azimu¬ 
thal  aperture,  20°;  the  two  eastward  windows 
("Tasman"  and  "Polynesian")  have  apertures  of 
8°  each. 

The  aperture  of  the  westward  window  is 
bounded  by  Brazil  and  West  Africa.  Sound 
transmitted  from  Heard  Island  has  unobstructed 
access  to  receivers  along  a  wide  swath  of  the 
northwest  Atlantic  Ocean  between  Nova  Scotia 
in  the  north  and  Virginia  in  the  south.  Acoustic 
signals  on  a  more  complicated  and  longer  path, 
from  Perth  to  Bermuda,  were  measured  during  a 
history-making  1960  experiment.  Interpretation 
of  these  measuremenLs  was  given  by  Schockley 
et  al.  (1982)  and  Munk  et  al.  (1988). 

The  main  eastward  window  is  azimuthally 
limited  to  the  "Tasman  Gap,"  an  opening 
between  Tasmania  and  North  Island  of  New 
Zealand.  It  may  be  further  limited  by  bathy¬ 
metric  barriers  created  by  the  Western  Polyne¬ 
sian  Islands  or  the  Hawaiian  Islands.  A  rela¬ 
tively  narrow  stretch  of  the  continental  slope 
along  the  North  American  Pacific  coast  off 
Washington  and  Oregon  seems  to  be  suited  for 
locating  receivers. 


The  o'her  eastward  window  is  bounded  by 
Wilkes  Land  (west  of  Addlie  Coast)  in 
Antarctica  and  South  Island  of  New  Zealand. 
Acoustic  paths  sweep  around  New  Zealand, 
pass  through  the  Polynesian  Archipelago,  and 
reach  the  continental  slope  off  California. 

While  the  geographical  "choke"  points 
that  bound  these  three  sets  of  paths  are  reason¬ 
ably  well  identified,  the  stretches  of  continental 
slopes  in  the  northwest  Atlantic  and  northeast 
Pacific,  where  we  expect  favorable  listening 
conditions,  can  be  adequately  identified  only  by 
modeling  the  refracted  acoustic  paths. 

These  are  some  of  the  longest  underwater 
propagation  paths  on  the  globe.  When  modeled 
with  realistic,  nonsmooth,  oceanic  features  of 
the  sound  speed  and  ocean  current  fields,  com¬ 
puted  ray  paths  over  such  long  distances  may 
exhibit  chaotic  behavior  similar  to  that 
described  by  Palmer  et  al.  (1988).  Chaos,  in  the 
analytical  and  computational  sense,  is 
suppressed  by  smoothness  in  the  representation 
of  the  ocean,  but  this  suppression  degrades  the 
realism  of  the  model.  Chaos  is  accentuated  by 
dimensionality  of  the  model,  in  that  3-D  trajec¬ 
tories  arc  much  more  prone  to  exhibit  chaos 
than  2-D  trajectories.  But  3-D  modeling 
enhances  the  realism  and  resolving  power  of  the 
model.  Hence  the  model  design  criteria  arc 
conflicting.  We  have  paid  attention  to  this 
important  topic  in  the  Discussion  (Sections), 
and  have  interspersed  throughout  the  text  con¬ 
siderations  that  will  facilitate  development  of 
"2-D-I-"  boundary  layer  ray  tracers  for  use  with 
nonsmooth  ocean  data.  However,  our  computed 
results  arc  obtained  with  a  3-D  Hamiltonian  ray 
tracer  in  a  moderately  smooth  ocean,  and  they 
do  not  exhibit  any  sign  of  exponential  sensi¬ 
tivity  to  initial  conditions  even  at  the  extreme 
range  of  18  Mm. 

In  the  work  reported  here,  our  objective 
was  to  explore  fundamental  propagation  condi¬ 
tions  along  the  "Tasman"  path  and  develop 
modeling  software  to  support  and  interpret  the 
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Heard  Island  feasibility  and  long-term  acoustic 
monitoring  experiments. 

A  SOFAR  acoustic  propagation  path  from 
Heard  Island  through  the  Tasman  Gap  to  the  US 
West  Coast  is  viable  (in  the  ray  approximation) 
if  the  SOFAR  axis  is  not  obstructed  by  bathy¬ 
metry.  We  do  not  address  bathymetry  effects  in 
this  report,  and  we  assumed  the  viability  of  a 
Tasman  path  from  the  results  of  Munk  and 


Forbes  (1989) — their  Fig.  6.  Our  ovm  analysis 
of  the  SOFAR  geometry  along  the  hypothetical 
Tasman  path  has  shown  their  Fig.  6  to  be  in 
error,  and  this  puts  the  viability  of  a  Tasman 
SOFAR  path  into  question.  Therefore,  the  com¬ 
puted  results  that  we  present  for  propagation 
along  the  Tasman  path  are  more  appropriate  for 
illustrating  modeling  techniques  than  for  draw¬ 
ing  geophysical  inferences. 


Figure  I .  Geodesics  bounding  the  westward  Atlantic  window  from  Heard  Island. 
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Se''!ion  Preview 

Section  2  is  the  main  body  of  this  report. 
It  addresses  geodesic  considerations  and  the 
Hamiltonian  formalism.  Section  3  describes  the 
models  used  to  represent  the  sound  speed  struc¬ 
ture  of  the  Indian/Pacific  Ocean  and  the  flow 
field  of  the  Antarctic  Circumpolar  Current.  Sec¬ 
tion  4  presents  general  modeling  results  and 
specific  eigenpatns  to  an  arbitrary  point  near 
Coos  Bay,  Oregon.  Section  5  is  a  discussion  of 
our  conclusions  and  touche.,  on  considerations 
for  the  design  of  new  modeling  software. 

2.  Earth  Coordinates  and  the  3-D  Hamil¬ 
tonian  Formalism 

Exp  \,iice  with  three-dimensional  under¬ 
water  acoustic  ray  tracing  over  ranges  oa  the 
scale  of  one  radius  of  the  earth,  or  larger,  is  very 
recent  and  not  generally  available  from  the 
literature.  When  reported,  the  results  are  typi¬ 
cally  based  on  substantial  simplifying  assump- 
dons. 

What  "Fully"  J-D  Means  to  Ray  Tracing 

The  results  shown  here  were  obiainrd 
from  fully  3-D  models,  and  it  is  important  to 
clarify  wnat  is  meant  by  that.  Ray  paths  that  are 
confined  to  a  (general)  cylindrical  surface  (with 
a  vertical  generatrix)  are  3-D,  but  they  are  not 
"fully  3-D"  because  they  are  confined  to  a  sur¬ 
face.  This  surface  is  "ruled";  i.e.,  it  is  particu¬ 
larly  simple  and  regular.  The  directrix  (the 
curve  through  which  the  generatrix  always 
passes)  of  such  a  cylinder  is  a  groat  circle,  a 
geodesic,  or  a  refracted  geodesic.  The  particu¬ 
lar  type  of  directrix  depends  on  the  sophistica¬ 
tion  of  tlie  given  model,  but  all  types  have  this 
in  comiiion:  at  all  "down-range"  points,  the  ray 
paths  have  the  same  azimuth  regardless  of  their 
launch  elevation  angle.  The  converse  is  also 
true.  V/hen  the  rays  launched  with  the  same 
azimuth,  but  different  launch  elevation  angles. 


exhibit  different  azimuths  at  an  arbitrary  down- 
range  point,  then  the  corresponding  ray  paths 
could  not  have  been  confined  to  a  mled  cylindr¬ 
ical  surface;  these  ra''  paths  are  fully  3-D. 

How  important  is  the  "fully  3-D"  ray  trac¬ 
ing?  i-jmch  (1990)  reports  that  the  out-of-planc 
refraction  effects  are  sn,all,  except  when  the  ray 
pathi  get  deflected  from  the  cylindrical  surface 
b;  reflection  (from  baih'^Tnetry)  or  by  strong  but 
spatially  uneven  refraction.  We  concur. 

Specifically,  the  propagation  geometry  of 
paths  that  oscillate  only  a  few  hundred  meters  in 
the  vertical  (about  the  sound  channel  axis)  can 
be  thought  of  as  effectively  confined  to  a  single 
cylindrical  surface.  A  horizontally  refracted 
geodesic  is  a  good  modeling  directrix  for  this 
cylinder.  Nonethelc:. ,  we  will  show  that  evtm 
for  such  paths  our  stringent  eigenray  tolerances 
cannot  be  met  without  capturing  the  azimuthal 
deviations  caused  by  out-of-plane  refraction. 

The  "fully  3-D"  modeling  that  we  actually 
used  may  be  much  more  important  for  rays  that 
oscillate  vertically  several  kilometers,  but  we 
have  not  modeled  such  rays. 

Global  and  Local  Coordinates 

In  underwater  acoustics,  there  is  a 
geometrical  advantage  in  that  the  ocean  is  a 
thin,  sharply  bounded  waveguide  with  a  thick¬ 
ness  less  than  1/1000  of  the  earth’s  radius.  This 
nearly  decouples  any  local  reference  system 
from  a  global  reference  system.  In  this  context 
"local"  is  logically  related  to  "vertical,"  and 
"global"  to  the  "horizontal"  extent  of  the  propa¬ 
gation  path. 

To  focus  attention  to  tbe  issues,  visualize  a 
"flat  earth"  modeling  (i.e.,  without  embedding 
curvature  terms  into  the  ray  equations)  of  under¬ 
water  propagation  on  ?  spherical  earth  to,  say, 
6000  km  range.  We  stress  the  existence  of  two 
types  of  errors.  The  first  would  be  the  inability 
of  the  modeled  ray  paths  to  curve  horizontally 
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along  a  great  circle  route — an  inability  caused 
by  the  intrinsic  lack  of  azimuthal  control.  As  a 
result,  the  modeled  sound  paths  would 
encounter  waters  of  substantially  different  pro¬ 
perties,  and  this  would  cause  travel-time  enors 
in  addition  to  the  errors  due  to  the  mis¬ 
represented  geometry. 

The  second  error  type  of  the  "flat  earth" 
model  is  due  to  the  earth’s  curvature  and  occurs 
in  the  vertical  section.  It  is  self-evident  there, 
but  one  should  keep  in  mind  that  the  computed 
underwater  sound  paths,  just  like  the  actual 
ones,  cannot  escape  the  oceanic  waveguide  no 
matter  how  bad  the  representation  of  the  propa¬ 
gation  geometry.  The  modeled  sound  will  reach 
the  targeted  di.stance  with  some,  but  not  catas¬ 
trophic,  error  in  arrival  time. 

The  "horizontal"  error  can  be  attributed  to 
an  inadequate  topology  of  the  Cartesian  system 
in  the  global  context  of  the  problem.  The  "verti¬ 
cal  .section"  error  reflects  accumulation  of  errors 
due  to  an  inadequate  metric  of  the  Cartesian 
system  in  the  local  context.  We  note  that 
mathematicians  would  call  both  deficiencies  of 
the  Canesian  sy.stem  "metric"  in  the  spherical 
context,  but  we  use  the  term  "topological"  to 
emphasize  a  difference  in  the  overall  shape  of 
the  propagation  paths.  In  the  horizontal  and  on 
a  large  scale  this  difference  would  be  quite  dis¬ 
cernible,  but  in  the  vertical  the  correct  and  the 
erroneous  paths  would  look  rca.sonably  alike  on 
a  large  scale,  as  they  zigzag  between  the  sea 
surface  and  the  ocean  bottom. 

How  to  "Flatten"  the  Earth's  Curvature 

Both  types  of  errors  are  as.sociated  with 
the  earth’s  curvature.  Actually,  only  the  second 
was  recognized  early;  it  was  dealt  with  using 
approximations  that  were  appropriate  given  the 
questions  asked  and  the  computational  limita¬ 
tions  at  the  time. 

For  ranges  up  to  one  earth  radius,  the 
second  (vertical)  effect  of  earth’s  curvature  can 
be  reasonably  accounted  for  by  introducing  a 


number  of  "earth  flattening"  devices.  The  best 
known  among  them  was  introduced  by  Booker 
and  Walkinshaw  in  1946.  It  is  a  modified  index 
of  refraction,  and  the  computing  is  still  done  in 
a  Cartesian  system.  The  validity  of  the  ensuing 
approximation  has  been  analyzed  by  Bre- 
khovskikh  (1960),  and  by  Budden  (1961;  1988) 
for  radio  waves  in  the  ionosphere. 

In  underwater  acoustics,  the  sound  speed 
profile  can  be  adjusted  to  compensate  for  the 
earth’s  curvature.  This  has  been  done  by  Wein¬ 
berg  (1981)  in  the  Generic  Sonar  Model  (GSM), 
closely  paralleling  the  modified  index  idea. 

The  practice  of  modifying  the  medium’s 
description  to  avoid  dealing  with  more  complex 
coordinates  is  probably  waning.  Brekhovskikh 
no  longer  refers  to  it  in  1980.  The  GSM  models 
a  layered  (range-independent)  medium,  and  this 
limits  the  distances  of  practical  interest. 

Beyond  the  Flat  Earth:  The  Perfect  Sphere 

.Most  long  distance  ray  tracing  near  and 
below  the  earth’s  surface  is  now  done  in  spheri¬ 
cal  coordinates.  This  removes  both  topological 
and  metric  errors  as  long  as  the  spherical 
representation  is  considered  to  be  sufficient  for 
the  problem  at  hand.  Among  the  models  that 
use  spherical  coordinates  in  underwater  acous¬ 
tics  arc  the  already  mentioned  HARPO  by  Jones 
ct  al.  (1986)  and  its  extension  (HARPO-MODl) 
by  Lynch  (1990).  An  overview  of  numerical 
ocean  acoustic  propagation  in  three  dimensions 
is  expected  soon  in  book  form,  from  Lee  (1990). 

The  ray  trace  equations  for  spherical  coor¬ 
dinates  arc  published  in  a  number  of  texts  and 
papers,  but  rarely  as  completely  and  computa¬ 
tionally  oriented  as  in  Jones  et  al.  (1986).  An 
accessible  account  of  ray  tracing  in  spherical 
coordinates,  as  used  in  seismic  work  a  decade 
ago,  is  given  by  Aki  and  Richards  (1980). 

Beyond  the  Sphere:  The  Ellipsoid 

The  use  of  a  3-D  spheroidal  coordinate 
system,  and  approximations  to  it,  must  be 
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motivated  by  need — the  need  being  to  eliminate 
the  error  due  to  the  horizontal  lay  of  acoustic 
Fermat  paths  over  very  long  distances  if  they 
are  modeled  as  great  circles.  The  latter  is,  of 
course,  dircctly  linked  to  ray  path  computations 
in  a  spherical  coordinate  system.  At  long  pro¬ 
pagation  ranges,  the  lay  of  Fermat  paths  is  sys¬ 
tematically  different  on  an  ellipsoid  (sphenoid) 
than  on  a  sphere.  Implications  of  this  fact  were 
thoroughly  analyzed  by  Munk  et  al.  (1988). 

The  longest  paths  for  underwater  sound  on 
this  planet  must  involve  the  southern,  near- 
Antarctic,  Indian  Ocean  no  matter  how  high  the 
latitude  of  a  receiver  or  transmitter  in  the  North 
Atlantic  or  North  Pacific.  The  intrinsic  ellip¬ 
soidal  topology  causes  the  true  Fermat  path  at 
such  ranges  to  always  project  "south"  of  the 
corresponding  great  circle  path.  This  argument 
is  based  on  differential  geometry,  but  remains 
valid  in  the  presence  of  arbitrary  lateral  gra¬ 
dients  of  the  sound  speed  field  along  the  path. 

An  oceanographic  consequence  of  this 
systematic  deviation  is  that  sound  will  traverse 
through  "colder"  (slower)  waters  than  it 
apparently  would  have  if  the  propagation  had 
been  modeled  in  spherical  coordinates.  Hence 
the  interest  of  Munk  et  al.  (1988),  as  well  as 
Munk  and  Forbes  (1989),  in  identifying  the  need 
for  modeling  propagation  as  close  as  topologi¬ 
cally  practical  to  the  geometry  of  the  planet. 

The  following  discussion  addresses  details 
of  the  geometrical  concepts  related  to  the  coor¬ 
dinate  .system.  The  salient  points  are 

•  The  earth’s  ellipsoid  offers  a  better  topology 
for  modeling  long  range  propagation  than 
does  a  sphere. 

•  The  ellipsoid  has  a  weak  eccentricity,  and  its 
metric  is  not  very  different  from  that  of  a 
sphere. 

•  The  ellipsoidal  coordinates  are  more  complex 
than  the  spherical  ones.  The  differential 
operators  have  more  terms  and  have  the 
eccentricity  implicitly  embedded  in  a  compli¬ 
cated  way.  If  they  are  evaluated  at  every 


integration  step,  they  may  become  a  computa¬ 
tional  monster,  but  if  omitted  the  coordinate 
system  becomes  an  approximation. 

•  The  ocean  data  no  longer  possess  the  cus¬ 
tomary  structure  if  they  are  to  be  strictly 
represented  in  the  ellipsoidal  coordinates. 
For  example,  data  given  at  a  constant  depth, 
but  different  position,  no  longer  have  the 
same  eUijjsoidal  "shell"  coordinate. 

•  The  maximum  ratio  of  vertical  ray  path  dis¬ 
placement  in  the  5-6  km  ocean  depth  to  the 
earth  radius  is  small  Oess  than  1  x  10"^),  and 
the  relative  change  in  the  spheroid’s  eccentri¬ 
city  over  this  displacement  is  equally  small. 
This  permits  a  local  representation  of  data 
(depth  can  be  used  instead  of  the  ellipsoidal 
shell  thickness),  and  the  inherently  local 
(differential)  computation  of  ray  position  can 
be  accurately  based  on  the  depth  coordinate, 
as  will  be  shown  in  the  discussion  on  the  radii 
of  curvature. 

The  Oblate  Spheroid 

The  earth’s  spheroid,  or  ellipsoid,  is 
mathematically  classified  an  oblate  spheroid. 
This  shape  is  an  ellipsoid  of  revolution,  revolv¬ 
ing  about  its  minor  axis.  It  possesses  less  sym¬ 
metry  than  a  sphere  and  thus  has  more  compli¬ 
cated  relations  between  its  coordinates. 
Figure  3(a)  (adapted  from  Baier,  1972)  gives  a 
good  3-D  picture  of  the  coordinate  surfaces  on 
the  spheroid.  The  surfaces  labeled  U3  =  const 
arc  the  "radial"  shells,  and  the  hyperboloidal 
mantle  surfaces  U2  =  const  are  "latitudes."  The 
coordinate  surfaces  U\  =  const  are  planes 
through  the  axis  of  rotation  and  are  defined  as 
"meridians." 

A  cross  section  through  the  rotational  axis 
of  such  an  ellipsoid  is  shown  in  Fig.  3(b).  The 
meridional  surfaces  are  no  longer  shown. 
Traces  of  the  U2  =  const  surfaces  (constant  lati¬ 
tudes)  are  hyperbolas,  and  they  must  be  orthog¬ 
onal  to  the  M3  =  const  shells  (though  they  are 
not  shown  as  such  in  this  textbook  sketch). 
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Figure  3.  (a)  Coordinate  surfaces  on  the  spheroid  [from  Fig.  I  in  "Acoustic  radiation  impedance 

of  caps  and  rings  on  oblate  spheroidal  baffles"  by  R.V.  Baier,  J.  Acoustic  So.  Am.,  5U5J.- 
1705-1716,  1972];  (b)  a  cross  section  through  the  rotational  axis  [from  Fig.  2.12  in 
Mathematical  Methods  for  Physicists  (2nd  ed.)  by  George  Arfken,  copyright  1970  by 
Academic  Press,  New  York]. 
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Details  about  these  coordinates  can  be 
found  in  Arfken  (1970)  (from  which  Fig.  3(b) 
was  drawn),  or  from  Morse  and  Feshbach 
(1953)  in  even  greater  detail.  Both  works  use  a 
left-handed  system  in  ordering  the  triplets 
(u3,M2.«i).  and  we  retain  it  (cyclically  per¬ 
muted  when  necessary). 

When  attention  is  confined  to  the  surface 
of  the  spheroid,  «3^,  then  the  geodesic  literature 
offers  a  wealth  of  detailed  and  specialized  infor¬ 
mation  (e.g.,  Bomford  (1980),  Zakatov  (1962), 
and  Jordan  et  al.  (1965).  The  relations  of  direct 
interest  for  long  range  acoustic  propagation  are 
not  too  readily  extracted  from  that  literature,  so 
they  will  be  quoted  or  derived  below  for  ease  of 
reference. 

Relations  Between  Spheroidal  Coordinates 

The  coordinates  of  an  oblate  spheroid  are 
usually  found  in  two  forms.  One  (Fig.  3) 
employs  the  left-handed  orthogonal  triad 
{uj,u2,u\},  for  which  one  can  express  the 
right-handed  Cartesian  triad  {x,y,z}  in  terms  of 
hyperbolic  functions  of  W3  and  trigonometric 
functions  of  M2  and  u\.  This  form  is  given  by 
Morse  and  Feshbach  and  by  Arfken. 

The  other,  completely  equivalent  form  is 
more  algebraic  (Madelung,  1957).  It  employs 
the  left-handed  orthogonal  triad  fK2,X2,'ki}.  It 
also  results  in  notationally  simpler  differential 
expressions.  This  is  helpful  while  inspecting 
derivations;  e.g.,  left-handedness  of  the  curvi¬ 
linear  coordinate  triads  can  be  tracked  by  apply¬ 
ing  a  minus  sign  to  the  Vx  operator. 

The  following  equivalences  hold  (cf. 
Fig.  3); 

X3  =  sinh  M3 

X2  =  sin  M2 

Xi  =  M J. 

Ranges  of  the  arguments  as  well  as  the 
interpretations  of  surfaces  generated  by  them 
are 


•  The  surfaces 

M3  =  const,  0<M3  <00 

are  oblate  spheroidal  shells  (so  that  M3  partly 
describes  a  "radial"  distance  to  the  shell’s  sur¬ 
face). 

•  The  surfaces 

M2  =  const,  -nil  <  U2  ^  Jc/2 

are  one-shelled  hyperboloids.  A  limiting  (and 
degenerate)  case  is  the  equatorial  plane, 
M2  =  0.  The  coordinate  M2  becomes,  by 
definition,  the  geodetic  latitude,  X- 

•  The  surfaces 

M  ]  =  const,  0  <  M 1  <  2jt 

are  half  planes  through  the  z  axis  (the  axis  of 
rotation),  and  u  1  is  identical  to  the  geodetic 
longitude,  (j). 

Relations  of  these  curvilinear  coordinates 
to  the  Cartesian  coordinates  are 

X  =  a<  COShM3  COSM2  COSMi 
=  a,  (1  (1  -  cos^i  ;  (1) 

y  =  a,  coshM3  COSM2  sinu  1 
=  a,  (1  (1  -  ^2^)'^^  sinXi  ;  (2) 

z  =  a,  sinhM3  sinM2 

=  a,  /\3  X2  .  (3) 

The  proportionality  factor  ot,  is  obtained 
by  elimination  of  coordinates,  using  the  follow¬ 
ing  geometrical  facts;  the  smaller  axis,  h,  of  the 
spheroid  satisfies  z^  =  for  X2  =  1,  and  the 
larger  axis,  a,  satisfies  =  a}  for  X2  =  0. 

Therefore,  a,  =  -  b^,  and  a,  has  a  geometri¬ 

cal  interpretation  as  the  radius  of  the  "focal  cir¬ 
cle"  of  the  oblate  spheroid  in  the  equatorial 
plane  (Fig.  3(a)). 

The  geocentric  radius  (distance  ftom  the 
origin  to  a  point  on  any  of  the  spheroidal  shells) 
is  important  for  the  metric  properties  of 
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approximate  computations  involving  the  earth’s 
ellipsoid;  it  is  given  by 

pg^=x^  +y^  +  z^ 

+  1)-  (4) 

Numerical  Geometry  of  the  Ellipsoid 

A  frequently  used  configuration  of  the 
earth’s  ellipsoid  is  parameterized  by  definitions 
of  the  WGS  84  (Department  of  Defense  World 
Geodetic  System  1984)  spheroid,  and  we  use  it 
here  in  our  text.  However,  most  of  the  com¬ 
puted  results  used  parameters  from  the  prede¬ 
cessor  WGS  72.  The  72  spheroid  is  negligibly 
different  from  the  84  (its  major  semiaxis  is 
shorter  by  2  m,  and  its  "flatness"  is  smaller  by 
3x10"*).  Two  parameters  are  sufficient  to 
define  the  geometry:  the  major  semiaxis  length 
a  =  6378.137  km  and  the  dimensionless  flatness 

/=  1/298.2572. 

From  a  and  f  one  obtains  the  minor  semiaxis 
b  =  and  the  "first"  or  "major"  eccentri¬ 
city  squared  e^  =  =  1/149.3790.  Note 

that  all  ellipsoids  in  use  so  far  (e.g.,  Bessel 
(1841),  Clarke  (1866),  International,  or  Hayford 
(1924))  have  their  ellipticity  described  within 
limits  of  the  approximations  /=1/3(X)  and 
e^  =  1/150. 

The  following  identities 

(l-/)2  s  l-c^  H  (6/a)2,  (5) 

in  conjunction  with  Eqs.  (1)  through  (4),  yield  a 
definition  of  the  earth’s  oblate  spheroid  terms: 

e^  (6a) 

(6b) 

uo  =  sinh"'  X3, .  (6c) 

where  the  subscript  0  indicates  the  parameter 

value  at  the  earth’s  surface.  Using  the  numeri¬ 
cal  parameters  of  WGS  84,  one  obtains 


ct,  =  521.85400842262  km 

(7a) 

X3,  =  12.18109320164 

(7b) 

uo  =  3.1947128 . 

(7c) 

The  extended  precision  of  values  in  Eq.  (7)  is 
carried  along  solely  for  checking  the  integrity  of 
the  approximate  computations;  their  accuracy 
depends  on  the  accuracy  of  the  ellipsoidal  axis 
a  (seven  significant  digits)  and  the  flattening  / 
(about  seven  significant  digits).  The  parameters 
for  our  earth’s  oblate  spheroid,  Eq.  (7),  could 
not  be  found  in  the  topical  literature,  probably 
because  geodesy  (including  the  physical  one) 
does  not  find  much  use  for  the  earth’s  geometri¬ 
cal  ellipsoid  at  any  distance  below  the  sea  sur¬ 
face.  However,  the  oblate  spheroidal  coordi¬ 
nates  are  useful  in  geophysical  contexts,  spurred 
by  the  use  of  earth-orbiting  satellites  (Vinti, 
1959). 

Approximations 

In  our  application,  the  fundamental  shape 
parameters  of  (7)  are  essential  for  high  reso¬ 
lution  checks  on  the  accuracy  of  approximate 
computations  when  ray  paths  dip  a  few  kilome¬ 
ters  below  the  sea  surface. 

A  few  qualitative  arguments  will  intro¬ 
duce  the  approximations  relevant  to  the  global 
oceanic  boundary  layer.  It  is  evident  from 
Fig.  3(b),  and  Eqs.  (1)  through  (3),  that  constant 
latitude  surfaces  on  the  spheroid  do  not  intersect 
the  rotational  axis,  unlike  the  constant  latitude 
surfaces  on  the  sphere  where  they  are  cones 
with  a  single  vertex  at  the  earth’s  center.  More 
importantly,  the  local  normal  to  any  X.3  =  const 
shell,  tangential  to  a  X.2  =  const  surface,  does 
intersect  the  rotational  axis  but  not  at  the  earth’s 
center. 

Which  Latitude? 

In  geodesic  terminology,  the  discussion 
above  gives  rise  to  two  latitudes,  the  geocentric 
latitude 
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V  =  tan"’  [z  /  (x^  +  , 

and  the  geodetic  latitude,  x.  introduced  previ¬ 
ously.  They  are  related  by 

tan\(i  =  (1  -e^)tanx, 


+  0.1 


^2^  +  ^3^ 
1-^2^ 


dXz^ 


+  aU^+h^)(^-h^)d\i^ .  (8) 


which  shows  that  x  >  V  everywhere  except  at 
the  poles  and  the  equator  where  they  have  equal 
values.  The  difference  between  them  is  largest 
at  latitude  45°  where  the  angle  is  11 '33".  In 
acoustic  propagation  modeling  this  difference 
can  be  interpreted  as  a  departure  from  the  "verti¬ 
cal"  that  affects  the  given  data  fields  (though  the 
concepts  of  vertical  and  normal  have  very  pre¬ 
cise  references  that  are  of  no  concern  here).  In 
that  sense  it  is  a  negligible  difference.  We  do 
not  distinguish  between  the  two  latitudes, 
though  we  do  follow  the  more  common  practice 
and  preferentially  use  the  geodetic  latitude,  X- 


Which  Radius? 


A  second  point  evident  in  Fig.  3(b)  is  that 
while  the  earth’s  surface  is  associated  with  the 
constant  value  of  the  spheroidal  coordinate  >,3,, 
which  would  conventionally  be  the  "radial" 
coordinate  in  a  spherical  analogue,  the  ellip¬ 
soidal  earth’s  surface  is  obviously  not  associ¬ 
ated  with  the  same  value  of  the  geocentric 
radius,  Eq.  (4),  at  arbitrarily  chosen  surface 
points.  In  the  pure  spherical  system,  all  surface 
points  have  the  same  geocentric  point-radius. 
The  neighborhood  of  each  point  has  a  unique 
radius  of  curvature  (identically  equal  to  the 
point-radius,  and  with  a  center  coinciding  with 
the  point-radius  origin)  that  is  independent  of 
direction  (azimuth)  along  the  surface. 


The  situation  is  more  complicated  in  the 
spheroidal  system.  Consider  a  general  ray  point 
at  a  given  depth  below  the  ocean  surface  along 
the  acoustic  propagation  path.  To  get  a  handle 
on  consistent  differential  changes,  one  derives 
the  total  differential  of  any  3-D  path: 


ds 


2 


^2^  +  X,3^ 
1  +X32 


d'K'i^ 


The  first  term  of  Eq.  (8)  represents  a  differential 
"depth  equation,"  as  can  be  seen  from  the  fol¬ 
lowing  argument. 


Consider  a  fixed  latitude  and  longitude. 
Then,  Eq.  (8)  shows  the  relation  between  a 
given  change  in  >13  and  the  corresponding  depth 
change  dz  (because  this  rfr  is  a  change  along  the 
normal  to  the  ellipsoid). 


dz^ 


^2^  +  ^3^ 
1  +  X32 


dX-i^  . 


(8a) 


The  interesting  fact  is  the  presence  of  the  lati¬ 
tude  factor  ^2  in  Eq.  (8a).  This  shows  that  a 
surface  of  constant  ^,3  cannot  be  at  constant 
depth  in  a  global  sense.  Two  limiting  cases  are 
easy  to  check  out:  at  the  pole  X2  =  1  and 

dz^  =  a]  dX-i^  ;  (8b) 

at  the  equator,  X2  =  0,  so 

X  ^ 

dz^  =  al  £  dX^^  ■  (8c) 

1  +  X3^ 


The  fundamental  parameters  derived  in 
Eq.  (7)  become  useful  here.  Our  approximation 
in  the  Hamiltonian  differential  equations  will 
assert  that  (8a)  and  (8b)  are  negligibly  different 
at  any  relevant  latitude  for  depth  changes  over 
the  topmost  few  kilometers  near  the  sea  surface. 
'This  is  because  X3^^  is  very  large  compared 
with  unity.  Thus  the  case  is  made  for  abandon¬ 
ing  a  strict  accounting  of  X3  in  favor  of  the 
layer-by-layer  differential  analysis  that  monitors 
depth  changes  of  the  ray  path  without  referenc¬ 
ing  the  point-origin  of  the  coordinates.  This  is 
essential  for  an  efficient  representation  of 
environmental  data  (sound  Sfjeed,  flow  veloci¬ 
ties.  bathymetric  depths,  etc.). 
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Conversely,  the  spheroidal  surface’s  cur¬ 
vature  radii  that  a  ray  reaches  at  any  point  in  its 
trajectory  deserve  careful  attention.  At  any 
point  on  the  surface  of  a  sphenoid,  or  on  the  sur¬ 
face  of  an  interior  shell,  there  is  an  infinity  of 
radii  of  curvature,  and  they  are  functions  of  lati¬ 
tude  and  the  azimuth  at  each  point.  Two  of 
these  radii  define  all  others.  They  are  the  prin¬ 
cipal  radii  of  curvature,  the  meridional  radius  r„ 
(curvature  in  the  meridional  plane),  and  the 
prime  vertical  radius  r^  (curvature  in  the  plane 
of  the  "prime  vertical";  this  plane  is  perpendicu¬ 
lar  to  the  meridional  plane  through  the  surface 
point). 

The  principal  radii  of  curvature  are  com¬ 
pactly  expressed  by  introducing  an  auxiliary 
function,  w,  of  e  and  latitude  X‘. 

=  1  -  e^sin^x  ■  (9) 

so  that 


and 


tion  that  is  depth  dependent.  Since  the  factor  w 
in  Eqs.  (9)-{l  1)  is  close  to  unity,  the  true  curva¬ 
ture  radii  "at  depth"  are  closely  approximated 
when  the  surface  curvature  radii  are  decreased 
by  the  local  depth,  z,  of  a  ray  point. 

'•m(z)  =  r„-z 

r„iz)  =  r^-z.  (12) 

The  curvature  approximation  error  is 

sufficiently  small  for  modeling  propagation  in 
the  ocean  (at  a  depth  of  5  km,  it  is  2  x  lO"’  at  the 
equator  and  3xl0~^  at  latitude  45°),  and  it  is 
probably  small  enough  even  for  modeling  pro¬ 
pagation  at  subacoustic  frequencies  (at  a  depth 
of  20  km,  it  is  3x10“®  at  the  equator  and 
1 X 10“^  at  latitude  45°). 

It  follows  from  (11)  and  (10),  using  (9), 
that  everywhere  on  the  spheroid 

r.>r„,.  (13) 

Furthermore,  the  radius  of  curvature,  in  the 
direction  of  any  azimuth  a*  is  determined  by  a 
relation  due  to  Euler: 


Equations  (9)  through  (11)  are  valid  only 
on  the  earth’s  spheroid  surface  if  one  considers 
a  and  e  fixed.  This  is  the  normal  practice,  but 
conceptually  one  could  let  a  and  e  vary  as  func¬ 
tions  of  depth.  As  one  looks  at  spheroidal  shell 
surfaces  that  are  below  the  sea  surface,  the 
equatorial  radius  a  decreases  and  the  first  eccen¬ 
tricity  e  increases  with  the  same  relative  rate. 
This  can  be  seen  by  taking  the  logarithmic 
derivative  of  (6a),  because  the  left-hand  side 
(the  scale  parameter  a,)  is  an  invariant  for  our 
planet.  Equation  (6)  can  be  "driven  in  reverse" 
starting  with  the  surface  parameters  in  Eq.  (7), 
and  using  Eq.  (8a)  to  determine  for  any 
prescribed  z  or  dz. 

Alternatively,  one  leaves  a  and  e  fixed  at 
their  surface  value  and  considers  an  approxima¬ 


cos^tt;,  sin^tt/, 

Mr^  = - + - .  (14) 

rm  '•y 

Taking  the  mean  value  integral  of  r^^  over  aU 
azimuths,  one  gets  the  earth’s  mean  radius  of 
curvature  r  at  a  given  geodesic  latitude: 

^  '■y  • 

Using  (13)  with  (14)  it  is  clear  that 

rm^ra^<  r, 

and  that 

r„<r<r^. 

Figure  4  shows  the  deviations  of  r^  and  of  r„ 
from  r  as  a  function  of  latitude.  The  deviations 
are  largest  at  the  equator  where  the  difference 
between  the  two  principal  radii  of  curvature 
exceeds  42  km.  One  sees  that  the  spheroid 
behaves  like  a  sphere  at  the  poles  (having  a 
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Figure  4.  Deviations  of  and  of  x^from  r  as  a  function  of  latitude. 
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single  measure  of  curvature,  albeit  of  larger 
radius — nearly  6400  km — than  the  earth’s  aver¬ 
age  radius  of  6371  km).  But  the  spheroid  does 
not  behave  like  a  sphere  at  the  equator.  In  the 
tropics  the  azimuthal  steering  on  the  sphere  and 
on  the  spheroid  will  be  significantly  different. 
Long  acoustic  propagation  paths  cannot  be 
"very"  long  unless  they  cross  the  equator,  and 
therefore  this  difference  matters. 

Azimuthal  Steering  on  the  Spheroid 

Equation  (14),  and  a  text  passage  in  the 
NOAA  document  for  HARPO  (Jones  et  al., 
1986,  pp.  115-116),  provide  an  instructive 
motivation  for  clarifying  the  "steering"  issue: 
"For  some  long-range  ray  calculations  in  the 
ocean,  a  spherical  model  of  the  earth  may  not  be 
accurate  enough.  Some  applications  may 
require  ocean  models  to  be  expressed  in  geo¬ 
detic  (c.g.,  spheroidal)  coordinates,  which 
would  be  transformed  to  spherical  coordinates 
for  iay  tracing.  However,  for  paths  of  a  few 
thousand  kilometers  and  less,  we  would  recom¬ 
mend  using  a  spherical  computational  coordi¬ 
nate  system  with  W(l)  set  to  the  local  radius  of 
curvature  of  the  geoid  in  the  propagation  direc¬ 
tion." 

The  term  "ocean  models"  in  this  passage 
refers  to  the  analytical  models  representing  the 
refractive  index  or  current  velocity  data,  and 
W(l)  is  the  only  admissible  value  of  the  radius 
for  the  computational  sphere  in  HARPO. 

We  have  fundamental  problems  with  parts 
of  this  recommendation.  It  is  clear  that  the  ray 
paths  cling  to  a  "curved  earth"  because  the 
"data"  (including  the  top  and  bottom  boundaries 
of  the  oceanic  waveguide)  are  given  in  a 
specific  "curved  space,"  not  because  the 
differential  equations  are  cast  in  a  particular 
curved  coordinate  system.  It  is  a  question  of 
adjusting  mathematics  to  the  data  vs  adjusting 
data  to  the  mathematics.  The  cited  passage 
recommends  a  variant  of  adjusting  mathematics 
to  the  data  for  paths  of  a  few  thousand 


kilometers  or  less,  and  adjusting  data  to  the 
mathematics  for  longer  paths.  There  are  strong 
reasons  for  providing  mathematics  that  adjusts 
to  the  data  in  all  conceivable  cases,  regardless 
of  the  scale  of  the  problem. 

Since  that  part  of  the  recommendation  that 
leaves  data  "as  is"  and  uses  a  local  radius  of  cur¬ 
vature  (i.e.,  adjust  the  mathematics)  is  also 
easier  to  implement  in  HARPO,  we  explored  it 
Equation  (14)  furnishes  the  "local  radius  of 
curvature ...  in  the  propagation  direction,"  but 
this  direction  changes  even  on  a  perfect  sphere 
(i.e.,  the  path  is  orthodromic  and  not  loxo- 
dromic),  so  one  value  of  W(l)  is  hardly  enough. 
With  just  one  representative  radius  of  curvature 
along  the  path,  the  results  were  really  not  distin¬ 
guishable  from  a  pilot  run  on  a  perfectly  spheri¬ 
cal  earth  of  nominal  radius.  This  is  understand¬ 
able  in  the  context:  both  runs  were  done  in 
spherical  coordinates,  except  the  two  radii  were 
somewhat  different. 

Next,  by  modifying  HARPO’s  code,  we 
used  the  local  azimuth  along  the  path  to  recom¬ 
pute  the  azimuthally  dependent  local  curvature 
radius  at  every  step.  The  results  were,  again, 
hardly  distinguishable  from  the  spherical  pilot 
run,  but  a  close  inspection  of  propagation  times 
revealed  small  differences  between  the  nms. 
The  topological  vs  metric  aspects  of  long  range 
acoustic  ray  computations  had  become  clearer: 

•  As  long  as  there  is  only  one  radius  built  into 
the  differential  equations  associated  with  the 
Hamiltonian,  the  azimuthal  steering  will  be 
topologically  the  same:  on  a  golf  ball,  on  the 
earth,  or  on  Jupiter,  it  will  be  a  great  circle. 
Two  distinct  radii  are  needed,  (10)  and  (11), 
to  achieve  a  steering  that  results  in  a  varia- 
tionally  (geodesic)  minimal  or  stationary  tra¬ 
jectory. 

•  The  use  of  a  local  curvature  radius  "in  the 
direction  of  propagation"  is  metrically 
correct.  If  one  possesses  an  independently 
prescribed  geodesic  path,  such  as  the 
"refracted  geodesics"  of  Munk  et  al.  (1988), 
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then  the  computation  in  the  vertical  does 
better  when  a  variable  radius  of  curvature  is 
used  instead  of  a  fixed  radius  (e.g.,  equal 
volume  or  equal  area,  etc.).  However,  as 
explained  earlier,  the  major  representation 
errors  in  long  range  acoustic  propagation 
have  a  topological  rather  than  a  metric  origin. 

From  here  on,  we  can  speed  up  the  pace  of 
presentation  by  referring  to  the  literature, 
mainly  Munk  et  al  (1988)  and  Jones  et  al. 
(1986),  without  reproducing  their  analytical 
treatment. 

Horizontally  Refracted  Geodesics  on  the 
Spheroid 

A  complete  treatment  of  unrefracted  (sur¬ 
face)  geodesics  and  refracted  geodesics  (at  the 
depth  of  the  sound  channel)  is  given  by  Munk  et 
al.  (1988).  We  have  already  discussed  their 
assumptions.  In  each  set  of  their  equations  (Eq. 
(1)  for  the  unrefracted,  Eq.  (2)  for  the  refracted 
case),  there  are  three  equations.  Two  relatively 
simple  ones  describe  changes  in  latitude  and 
longitude,  and  the  more  complex  one  describes 
change  of  the  azimuth.  In  the  refracted  case, 
when  the  path  length  is  used  instead  of  the  pro¬ 
pagation  time  as  the  independent  argument,  the 
sound  speed  and  its  derivatives  occur  only  in 
the  azimuthal  equation.  This  reinforces  the 
notion  that  this  equation  is  the  (horizontal) 
"steering"  equation  for  the  propagation  path. 

The  following  topology  with  respect  to 
curvature  radii  is  relevant.  The  meridional 
radius  occurs  in  the  latitude  equation  (plausible 
and  easy  to  visualize),  the  prime  vertical  curva¬ 
ture  radius  occurs  in  the  longitude  equation,  and 
both  have  a  place  in  the  azimuthal  equation.  If 
the  sound  speed  were  constant  there  would  be 
no  difference  between  the  refracted  and 
unrcfracted  equations,  and  if  both  principal  cur¬ 
vature  radii  become  the  same  (as  they  nearly 
would  in  the  Arctic  Ocean)  then  everything 
reverts  to  the  spherical  case. 


Antipodal  Acoustic  Paths  on  the  Spheroid 

Munk  et  al.  (1988)  discuss  the  longest 
experimental  sound  propagation  path,  from 
Perth  to  Bermuda,  that  has  been  documented. 
The  range  was  19,820  km,  only  about  185  km 
short  of  antipodal.  They  are  also  interested  in 
true  antipodal  conditions  because  such 
geometries  may  be  realizable  in  the  future,  but 
this  is  not  a  primary  object  of  their  paper.  They 
say  that  the  "great  circle  geometry  fails  catas¬ 
trophically  for  near-antipodal  ranges,"  and  that 
"for  exactly  antipodal  transmissions  (geodesic 
180°)  the  geodesic  goes  through  the  pole,  as 
does  one  of  the  infinite  number  of  great  circle 
routes"  (their  Fig.  5).  These  are  correct  and 
important  observations,  particularly  relevant  in 
view  of  their  concluding  remark  that  this  land¬ 
mark  experiment  may  be  worth  repeating  given 
our  present  understanding  of  the  issues 
involved. 

We  are  specifically  interested  in  the  forth¬ 
coming  Heard  Island  feasibility  experiment. 
Heard  Island  is  closer  to  Africa  than  it  is  to  Aus¬ 
tralia,  and  the  acoustic  propagation  geometries 
from  Heard  Island  are  not  antipodal.  Heard 
Island’s  antipode  is  on  land  in  western  Canada, 
as  can  be  seen  in  Figs.  1  and  2.  However,  it  is 
clear  that  with  a  ship  under  way,  say,  from  Perth 
or  Freemantle  to  Heard  Island,  one  reaches 
potential  oceanic  transmission  points  that  have 
antipodal  geometry  relative  to  some  feasible 
receiver  locations  in  the  North  Atlantic. 

Two  questions  regarding  potential  under¬ 
water  acoustic  propagation  on  antipodal  paths 
arise  naturally; 

•  Are  antipodes  acoustically  "reachable"? 

•  If  reachable,  is  there  a  likelihood  of  enei]gy 
focusing  at  the  antipode?  If  not,  what  sort  of 
energy  focusing  can  be  expected  in  some 
neighborhood  of  the  antipode? 

Munk  et  al.  (1988)  show  that  horizontal 
refraction  significantly  modifies  the  actual 
sound  paths  relative  to  those  defined  by  spheri¬ 
cal  or  spheroidal  geodesics.  Such  modification 
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is  one  of  degree  rather  than  being  related  to  the 
underlying  fundamentals. 

TTius  a  conceptual  dilemma  exists: 

•  In  the  spherical  approximation  the  antipode  is 
reachable,  albeit  in  a  singular  way,  by  many 
great  circles  passing  through  a  feasible  oce¬ 
anic  window.  This  window  has  about  the 
same  azimuthal  width  as  the  corresponding 
window  in  the  spheroidal  approximation. 
Focusing  can  be  strong  (they  estimate  +20  dB 
for  the  Perth-Bermuda  geometry,  1 85  km 
short  of  the  antipode),  and  becomes  theoreti¬ 
cally  unlimited  at  the  antipode. 

•  Alternatively,  in  the  spheroidal  approxima¬ 
tion,  the  antipode  is  not  reachable  because  the 
unique  geodesic  to  the  antipode  must  go 
through  a  pole,  and  is  definitely  outside  the 
feasible  oceanic  window. 

We  investigate  these  questions  in  a  short 
aside  that  fits  the  geodesic  theme  of  this  section. 
Figures  5  and  6  describe  the  findings  for  a 
viable  propagation  geometry.  To  emulate  a  pro¬ 
pagation  window  of  similar  azimuthal  width 
(about  20°)  to  the  "Atlantic"  window  from 
Heard  Island  (Fig.  1),  we  place  a  hypothetical 
source  (Fig.  5)  at  the  longitude  of  Perth  but  at  a 
higher  latitude  (40°S).  The  horizontal  center  of 
the  radiated  fan  passes  through  Heard  Island, 
and  a  hypothetical  antipodal  receiver  is  in  the 
North  Atlantic  about  475  n.mi.  NNW  of  Ber¬ 
muda. 

Figure  6  shows  the  computed  pattern  in 
the  geographic  neighborhood  of  the  antipode. 
The  horizontal  ray  fan  focuses  on  an  area  about 
20  km  away  from  the  antipode,  which  is  not 
acoustically  reachable  through  the  ocean.  The 
focusing  is  not  to  a  point  but  is  spread  along  a 
caustic,  part  of  which  we  label  "audible"  in  the 
belief  that  the  pattern,  if  not  its  exact  location, 
will  persist  even  in  the  presence  of  lateral 
refraction. 

This  dilemma  was  clarified  by  Helmert 
(1880),  who  investigated  the  uniqueness  of  geo¬ 
desics  on  the  spheroid.  The  prototype  of  a 


surface  with  nonunique  geodesics  is  a  cylinder 
of  revolution.  Take  any  pair  of  distinct  points 
not  on  the  same  meridian,  the  eigen-geodesics 
between  these  two  points  are  helices  of  1, 2, ... , 
n,  and  up  to  an  infinite  number  of  windings. 
They  are  geodesics  because  each  has  a  shorter 
(stationary)  path  length  relative  to  any  of  its 
neighboring  paths.  The  helix  with  a  single 
winding  has  the  absolute  shortest  path  length;  it 
is  the  Fermat  geodesic.  The  spheroid  is  topo¬ 
logically  between  the  sphere  and  the  cylinder, 
but  our  spheroid  is  very  close  to  a  sphere.  Hel¬ 
mert  has  done  remarkably  detailed  computa¬ 
tions  of  the  geodesic  geometry  in  the  neighbor¬ 
hood  of  the  antipode  and  has  identified  the  mul¬ 
tipath  region  that  constitutes  our  present  con¬ 
text. 

The  caustics  that  we  label  "audible"  and 
"inaudible"  in  Fig.  6  constimte  one  leg  of  an 
astroid,  a  star-like  envelope  of  geodesics  cen¬ 
tered  on  the  antipode.  The  term  "inaudible" 
should  be  understood  in  context:  each  of  the 
four  legs  of  the  astroid  is  a  horizontal  caustic 
(i.e.,  an  envelope  curve  to  a  family  of  geo¬ 
desics),  and  only  part  of  one  of  them  will  be 
"audible."  The  family  of  geodesics  emanating 
from  the  Indian  Ocean  transmit  point  is  unique 
only  outside  this  astroid. 

The  major  axes  of  the  astroid  are  aligned 
N-S  and  E-W.  The  size  of  the  astroid  depends 
on  the  eccentricity  of  the  spheroid  and  on  the 
latitude  of  the  transmit  point.  In  particular,  the 
size  of  the  N-S  axis  is  the  greater  of  the  two, 
and  is  proportional  to  the  cosine  of  the  transmit 
latitude.  The  E-W  axis  size  is  proportional  to 
the  cosine  squared  of  the  transmit  latitude. 
Consequently,  high  latitude  transmit  points  are 
associated  with  a  small  region  of  geodesic  mul- 
lipathing  at  the  antipode.  The  size  of  the  astroid 
would  shrink  to  zero  if  a  geographic  pole  were  a 
feasible  transmit  point. 

Moreover,  the  relative  focusing  (+20  dB  at 
Bermuda)  estimated  by  Munk  et  al.  using  a 
spherical  argument  is  basically  correct  even  in 
the  spheroidal  case,  as  evidenced  by  Fig.  5. 
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Figure  5.  Geodesics  from  a  location  south  of  Perth.  The  fan  of  geodesics  with  azimuths  between 
220°  and  240°  approaches  the  Atlantic  antipode.  The  geodesic  with  an  azimuth  of  180° 
connects  the  antipodes  exactly  (181°  is  computable). 
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Figures  5  and  6  arc  generated  from  the  same 
high  precision  geodesic  computations,  but  the 
ray  fan  in  Fig.  5  is  very  similar  to  a  ray  fan  on  a 
spherical  earth.  Its  gross  structure  does  not 
even  hint  at  the  details  that  are  revealed  by 
Fig.  6. 

Hamiltonian 

We  chose  HARPO  as  a  preliminary  tool 
for  3-D  modeling  of  very  long  range  acoustic 
propagation.  It  was  e'ear,  though,  that  new 
software  would  eventually  be  needed  for  easier 
absorption  of  data  updates  on  a  global  scale. 
Then  the  ray  equations  could  be  carried  over  as 
Hamiltonian  equations  copied  from  HARPO,  or 
they  could  be  written  using  the  eikonal  equation 
as  a  base.  TTiis  accounts  for  our  interest  in 
understanding  the  Hamiltonian  formalism  in  a 
wider  setting. 

A  few  comments  about  the  Hamiltonian 
appioach,  in  parallel  with  those  given  by  Jones 
et  al.  (1986),  may  help  explain  why  numerical 
codes  that  arc  not  based  on  the  Hamiltonian 
may  serve  specific  objectives  rather  well. 

•  Not  only  arc  the  ray  paths  computed  from  the 
Hamilt'inian  equations  the  same  as  those 
computed  from  dilTcrential  equations  that  are 
derived  from  the  eikonal,  but  the  Hamiltonian 
iLsclf  is  derivable  from  the  eikonal  (Whitham, 
1961)  This  is  true  for  a  broad  formulation  of 
the  di.spcrsion  relation  (e.g.,  including  a  mov¬ 
ing  medium). 

•  It  is  also  true  that  the  Hamiltonian  is  the  more 
original  and  fundamental  construct  if  onv, 
thinks  of  the  eikonal  as  a  related  but  different 
concept.  This  relation  was  the  subject  of  a 
little-known  polemic  exchange  in  1937 
regarding  geometrical  optics,  but  applicable 
to  geometrical  acou.stics,  between  J.  L.  Synge 
and  M.  Herzberger.  It  ended  on  the  notion 
that  Bruns  (the  formulator  of  the  eikonal) 
rediscovered  in  1895  the  ray  propagation 
Hamiltonian  tliat  was  conceived  in  1832,  but 


brought  to  light  only  in  1931  (Conway  and 
Syng-). 

•  The  Hamiltonian  approach  to  geometrical 
acoustics  entered  the  literature  in  the  1950’s 
(Keller,  1954).  This  was  nearly  synoptic  with 
actual  applications  of  the  Hamiltonian  ray 
equations  in  radio-wave  propagation  by  J. 
Haselgrove  (1954).  Her  methodology  stimu¬ 
lated  much  of  the  subsequent  woilc,  bridging 
applications  in  several  media  of  propagation. 

•  The  literah  .re  cited,  and  the  individual  contri¬ 
butions  of  Jones  et  al.,  provide  a  good  exam¬ 
ple  of  some  activities  in  the  1960’s  and 
1970’s  leading  to  HARPO.  An  overview  by 
Ostashev  (1985)  puts  the  Hamiltonian  and 
eikonally  based  ray  formalisms  into  the  con¬ 
text  of  the  more  general  acoustic  propagation 
theory  for  an  inhomogeneous  moving 
medium. 

•  One  may  conclude  on  balance  that  the  Hamil¬ 
tonian  approach  offers  an  elegant,  but  not 
always  simple,  path  to  ray  trfcing.  This  is 
well  articulated  by  Lighthill  (1979)  who 
relegates  the  Hamiltonian  equations  to  what 
he  calls  a  "brief  parenthesis  on  parallel  with 
other  fields  of  study"  in  his  extensive  sections 
on  the  general  theory  of  ray  tracing  and  ray 
tracing  in  a  wind. 

•  1  he  primary  appeal  of  the  Hamiltonian  is  that 
it  is  independent  of  the  coordinates  used.  The 
converse  is  also  true;  the  Hamiltonian 
differential  equations  in  a  given  coordinate 
system  do  not  change  wh'^n  the  Hamiltonian 
changes  because,  say,  a  different  set  of  propa¬ 
gation  conditions  is  to  be  modeled.  We  rely 
heavily  on  this  in  modifying  HARPO’s  code 
for  use  in  the  ocean  on  the  ellipsoidal  earth. 
However,  the  ray  equations  derived  from  the 
eikonal  also  have  a  vector  form  (Pierce,  1981) 
that  is  intrinsically  independent  of  the  compu¬ 
tational  coordinate  system. 

There  are  six  ordinary,  first  order,  non¬ 
linear  Hamiltonian  equations  that  define  the  ray 


-18- 


APL-UWTR8929 


trajectory  in  a  three-dimensional  space.  Three 
of  them  define  the  position  coordinates  of  the 
ray  point.  Since  we  start  out  with  a  spherical 
coordinate  system,  all  terminology  refers  to 
such  a  system.  The  computed  ray-pomt  position 
is  the  (gcoccntrical)  radius  r,  the  colatitude  0, 
and  the  longitude  <().  The  other  three  equations 
define  the  local  components  of  the  wavenumber 
vector  1?  that  is  tangent  to  the  ray  path.  These 
six  equations  arc  integrated  subject  to  six  initial 
conditions;  the  three  coordinates  of  the 
transmitter,  and  the  three  components  of  the  ini¬ 
tial  direction  (initial  value  of  the  wavenumber 
vector). 

In  adapting  the  given  set  of  Hamiltonian 
equations  to  propagation  in  the  thin  ellipsoidal 
planetary  waveguide,  one  recognizes  immedi¬ 
ately  the  asymptotic  and  the  actual  structure; 

•  The  "horizontal"  components  of  the 
wavenumber  vector  (Kq.Ko)  must  implement 
the  az.imuthal  steering  that  corresponds  to  a 
"refracted  geodesic"  of  Munk  ct  al.  (1988). 
This  is  fed  back  to  the  latitude  and  longitude 
cqi  'on,  ensuring  a  correct  overall  global 
trajcL..iiy.  Four  of  the  six  equations  arc 
involved  in  this  process,  and  this  foursome  is 
strongly  coupled. 

•  The  other  two  equations  arc  the  radial  posi¬ 
tion  and  the  radial  component  of  the 
wavenumber.  One  more  equation,  a  horizon¬ 
tal  component  for  the  wavenumber,  is  needed 
to  define  the  ray  path  in  the  "vertical"  cylindr¬ 
ical  surface.  But  two  more  arc  available,  the 
aforementioned  equations  for  (Kq.k^).  They 
too  compri.se  a  foursome.  Again,  this  four¬ 
some  is  strongly  coupled. 

•  Thus  the  six  Hamiltonian  equations  subtend, 
in  a  topological  scn.se,  two  overlapping  four¬ 
somes  of  equations  with  strong  coupling 
within  each  foursome  and  weak  coupling 
across  them. 

•  The  weak  coupling  across  these  two  equation 
sets  is  implemented  in  the  numerical  code, 
and  it  ensures  that  the  ray  paths  arc  ful'y  3-D, 


i.e.,  not  conhned  to  an  idealized  cylinder  to 
which  they  would  be  bound  in  the  asymptotic 
case. 

The  ocean  data  (sound  speed  field  and 
current  velocity  field)  remain  defined  in  spheri¬ 
cal  coordinates  (latitude,  longitude,  and  depth). 
Depth  is  defined  as  the  difference  between  two 
geocentric  radii,  the  radius  to  the  sea  surface 
and  the  radius  to  the  point  considered.  For 
HARPO,  the  ocean  data  are  appropriatedly 
termed  "ocean  models"  because  they  must  be 
given  analytically  rather  than  pointwise. 

The  earth’s  ellipticity  is  introduced  into 
the  Hamiltonian  equations  by  distinguishing  the 
two  principal  radii  of  curvature,  r„  from  (10) 
and  Tj,  from  (11),  in  the  individual  terms  of  the 
equations  described  below.  The  depth  of  a  ray 
point  is  based  on  the  geocentric  radius  r  = 
from  (4).  Equation  (12)  may  then  be  substituted 
for  (10)  and  (11). 

A  tcrm-by-icrm  validation  of  the  modified 
equations  is  accomplished  by  requiring  that  they 
reduce  to  the  surface  geodesic  equations  when 
all  lateral  gradients  of  the  refractive  index  are 
.set  to  zero,  and  when  an  artificial  vertical  sound 
channel  is  established  with  its  axis  at  the  sur¬ 
face.  A  numerical  validation  was  performed  by 
running  this  reference  case  through  the  modified 
code. 

The  spherical  Hamiltonian  equations  that 
we  adapted  to  the  earth’s  ellipsoid  are  those 
residing  in  HARPO ’s  subroutine  HAMLTN. 
The  original  equations  are  (6.10)  through  (6.15) 
in  the  HARPO  document  by  Jones  et  al. 

Let  //  be  the  Hamiltonian,  t  the  propaga¬ 
tion  time,  Co  a  reference  constant  sound  speed, 
c  the  local  sound  speed,  p  =  co^the  independent 
variable  of  integration,  ©  a  fixed  wave  fre¬ 
quency,  and  F  =  co^W/^w  a  fixed  scale  factor. 
The  gcoccntrical  radius  to  the  ray  point  is  r,  the 
colatitudc  is  0,  and  the  longitude  is  Notice 
that  p  is  a  "scaled  time"  with  numerical  value 
roughly  equal  to  the  path  length.  The  exact 
relation  between  p,  t,  and  the  actual  path  length 
.V  is 


-19- 


APL-UWTR8929 
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The  six  new  Hamiltonian  equations  are 
dp  F  9k, 


M.  ^  1  dH 
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dp  r^sin0 


-  K<,r„cos0 


19//  .  a  dr 


dQ 


dp 


(15) 

(16) 
(17) 


(18) 


(19) 


(20) 


The  Hamiltonian  H  is  used  directly  as 
specified  by  Jones  et  al.  in  dispersion  relations 
designed  to  cover  specific  propagation  condi¬ 
tions.  Two  important  ones  are 

(i)  Propagation  through  an  inhomogeneous 
stationary  medium.  This  propagation  is 
isotropic;  the  ray  paths  are  perpendicular 
to  the  wavefronts. 

(ii)  Propagation  through  an  inhomogeneous 
moving  medium.  This  propagation  is 
anisotropic;  the  ray  paths  are  no  longer 
perpendicular  to  the  wavefronts. 


Ray  tracing  results  for  the  Heard  Island  to 
Oregon  Coast  paths,  utilizing  (15)  through  (20) 
and  both  isotropic  and  anisotropic  media,  are 
presented  in  Section  4. 

One  problem  related  to  data  representation 
remains  when  using  (15)  through  (20).  These 
are,  as  stated,  geocentrical  equations  on  the 
ellipsoid,  but  they  are  not  oblate  spheroidal 
equations.  For  them,  even  a  level  sea  surface, 
or  a  flat  bottom,  has  a  radial  position  that  varies 
with  the  horizontal  coordinates. 

Throughout  this  subsection  on  the  Hamil¬ 
tonian,  we  drew  attention  to  alternative 
viewpoints  that  are  relevant  to  ray  tracing.  A 
major  alternative  that  loomed  in  the  background 
is  the  abandonment  of  a  sharp  distinction 
between  2-D  and  3-D  in  the  thin  global  boun¬ 
dary  layer.  This  involves  a  "2-D-t-"  boundary 
layer  coordinate  system,  and  is  beyond  the 
scope  of  this  text.  However,  the  ray  tracing 
equations  will  be  looked  at  in  this  emerging 
context.  Therefore,  before  reporting  some  of 
the  numerical  results  obtained  using  a  3-D 
Hamiltonian  ray  tracer,  we  draw  attention  to 
Appendix  A.  In  it  we  illustrate  a  2-D  set  of  ray 
equations  in  two  coordinate  systems,  show  their 
relation  to  the  Hamiltonian  equations,  and  point 
to  the  merit  of  having  specialized  equations  that 
affect  possible  cancellations  analytically  rather 
than  numerically.  The  illustrated  equations  are 
those  that  Bold  and  Birdsall  (1986)  discuss  as 
the  "third  algorithm  (angular  forai),"  and  that 
are  used  in  a  number  of  interesting  applications 
(e.g.,  Baxter  and  Orr  (1982)  and  Eliseevnin 
(1965). 


3.  Indian/Pacific  Sound  Speed  Model  and 
the  Antarctic  Circumpolar  Current 
Model 

The  computer  code  in  HARPO  requires 
that  the  3-D  field  of  the  sound  speed  and  of  the 
current  velocity  be  continuous,  and  that  they 
have  continuous  partial  derivatives  in  the  three 
spatial  directions  (Jones  et  al.,  1986).  For 
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vector  fields,  each  component  must  satisfy  this 
requirement.  The  menu-driven  HARPO  further 
requires  that  such  fields  obey  a  particular  for¬ 
malism.  This  formalism  is  generally  termed 
OCEAN  MODELS,  and  it  includes  additional 
species  of  models,  for  the  ocean  bottom,  absorp¬ 
tion,  etc.  In  any  HARPO  run  there  may  be  no 
more  than  one  model  of  each  species,  but  some 
fields  allow  a  subspecies  of  the  general  model. 
The  general  model  for  a  field  is  termed  a  "back¬ 
ground"  model,  and  the  subspecies  is  termed  the 
"perturbation"  model  for  that  field.  This 
arrangement  permits  the  user  to  break  up  the 
complexity  of  modeling  a  field,  and  to  investi¬ 
gate  sensitivity  questions  by  switching  the  per¬ 
turbation  model  on  and  off. 

The  sound  speed  field  for  modeling  propa¬ 
gation  paths  in  the  Tasman  window  (Fig.  2)  had 
to  be  defined  over  ranges  that  were  so  long 
(18  Mm)  that  the  existence  of  a  single  perturba¬ 
tion  submodel  was  not  sufficient  to  describe  all 
significant  departures  from  the  background. 
Instead,  the  "background"  field  itself  had  to 
include  practically  all  the  complexity  that  was 
required.  Figure  7  shows  an  image  of  the  final 
analytical  sound  speed  field  as  a  function  of 
depth  (0  to  5  km)  and  latitude  (70°S  to  55°N), 
with  the  sound  speed  color  coded  between 
1450  m/s  and  1550  m/s.  The  rationale  for  this 
model,  and  its  limitations,  will  be  discussed 
after  pre.scnting  its  analytical  configuration. 

Sound  Speed  Model 

The  Indian/Pacific  Ocean  3-D  sound  speed 
model  is  valid  in  the  neighborhood  of  propaga¬ 
tion  paths  from  Heard  Island  through  the  Tas¬ 
man  window.  It  is  constructed  as  follows: 

•  In  the  vertical  coordinate  the  "background 
ocean  model"  satisfies  the  equations  of  a  con¬ 
tinuum  of  canonical  sound  speed  profiles 
(Munk,  1974)  and  is  referred  to  simply  as 
Munk’s  SSP  (sound  speed  profile).  This 
profile  was  primarily  meant  to  be  applicable 
in  tile  vicinity  of  the  sound  channel  axis,  but 


in  our  experience  it  can  match  archived  data 
over  a  significantly  wider  depth  range  at 
nearly  every  latitude  of  interest. 

The  functional  form  for  the  sound  speed 
c(z)  as  function  of  depth  z  (positive  down)  is 
defined  in  the  Munk  SSP  as 

c(z)  =  co[l -He(Y-i-e^-l)].  (21) 

where 

7(2)  =  (2-zo)- 

He 

This  profile  is  completely  determined  by  four 
parameters;  zq,  depth  of  the  sound  channel 
axis;  cq,  sound  speed  at  axis  depth;  e,  a  dimen¬ 
sionless  perturbation  or  shape  parameter,  and 
He,  a  characteristic  scale  depth. 

A  set  of  the  fitted  profiles,  offset  by  20  m/s 
each  and  indexed  by  latitude,  is  shown  in  Fig.  8. 
Table  I  contains  some  of  the  numerical  data  that 
were  used  to  estimate  the  four  parameters,  and, 
in  the  last  two  columns,  the  fitted  sound  speed  at 
the  surface  and  at  a  depth  of  4400  dbars.  The 
archived  data  used  to  support  this  model  were 
not  selected  by  "objective  means"  (least  square 
fits,  or  equivalents)  but  were  chosen  after  a  sen¬ 
sitivity  analysis  of  Eq.  (21)  that  is  outlined  next. 

Data  from  five  sources  were  used  in  this 
analysis.  The  conditions  at  the  sound  channel 
axis  were  obtained  from  the  appropriate  graphs 
(Figs.  3(a)  and  (b)  of  Munk  and  Forbes,  1989). 
For  the  rest  of  the  water  column  along  the 
acoustic  path,  the  sound  speed  was  computed 
from  oceanographic  data  extracted  from  three 
atlases  (Gordon  et  al.,  1982;  Wyrtki,  1988; 
Craig  et  al.,  1981).  Also,  tables  and  graphs  of 
Podeszwa  (1976)  were  used  to  estimate  the 
sound  speed  data  in  the  North  Pacific. 

A  few  interesting  facts,  not  described  in 
Munk’s  original  paper,  can  be  deduced  from 
Munk’s  canonical  profile.  They  show  the  versa¬ 
tility  of  Eq.  (21)  and  the  relevance  of  the  under¬ 
lying  physics: 
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•  The  parameter  pair  (e,  H^)  is  independent  of 

the  parameter  pair  (zq,  Cq)  because  3c/0e  0 

and  dddHc  — » 0  as  z  — »  zq.  Therefore,  the 
parameters  defining  conditions  at  the  sound 
channel  axis  can  be  chosen  independently 
from  those  that  define  the  upper  and  lower 
branch  of  the  SSP.  These  relations  are  illus¬ 
trated  in  Fig.  9.  Note  that  (zq.  Cq)  determines 
a  single  point,  while  (e.  He)  controls  the  inter¬ 
val  behavior  of  c(z)  away  from  the  sound 
axis. 

•  When  (zq.cq)  is  chosen  at  different  latitudes 
along  the  Tasman  window  propagation  path 
and  (e,He)  is  estimated  from  the  near  surface 
and  the  deep  water  sound  speeds,  then  He  and 
e  show  correlation  in  the  form  of  an  empirical 
affine  transformation; 

Fc[km]  =  aie  +  Oo  ,  (22) 

where  Oj  =115  and  Oq  =  0.3.  The  estimates 
ie,Hc)  from  Table  I  are  displayed  in  Fig.  10, 
and  it  is  seen  that  He  "tracks"  e  over  a  wide 
interval  of  latitudes.  The  scatter  diagram 
(e.,Hc)  on  which  Eq.  (22)  is  based  is  shown  in 
the  upper  right  comer  of  Fig.  10.  But  the 
apparent  relationship  implied  by  (22)  must  be 
interpreted  guardedly.  If  nature  behaved  this 
regularly  it  would  be  disappointing,  as  well  as 
unlikely. 

What  is  seen  in  Fig.  10  can  be  qualita¬ 
tively  explained  by  noting  another  asymptotic 
behavior  of  (21):  when  e  — » 0,  then  c  — »  Cq  for 
any  He  and  for  all  z — hence  the  term  "perturba¬ 
tion  parameter"  for  e.  This  parameter  attains  its 
largest  values  where  the  sound  channel  axis  is 
deep  (between  SO^S  and  40°S),  to  account  for 
the  sound  speed  contrast  between  surface  and 
deep  water  relative  to  the  water  at  the  channel 
axis.  TTie  scale  depth  He  must  also  be  large 
there  to  prevent  the  sound  speed  contrast  from 
becoming  too  large.  Conversely,  e  will  be 
small,  and  He  will  be  small  when  the  surface 
water  is  warm  (high  sound  speed),  but  the  chan¬ 
nel  axis  moves  closer  to  the  surface,  as  in  the 
low  latitudes  of  the  North  Pacific. 


This  explanation  shows  why  a  negative 
correlation  between  e  and  //c  as  e  or  evolves 
along  a  hypothetical  propagation  path,  could 
never  fit  the  SSP  of  any  deep  ocean;  a  small  e 
and  large  He  would  estimate  an  SSP  that  is  too 
uniform  top  to  bottom,  while  a  large  e  and  small 
He  would  imply  sound  speeds  that  are  much  too 
high  at  the  surface,  or  at  the  bottom,  or  both. 
Thus,  choices  of  e  and  He  that  smoothly  track 
the  path  parameter,  but  randomly  map  about  a 
positive  correlation  line,  can  help  generate 
canonical  profiles  that  effectively  model  the 
vertical  distribution  of  the  sound  speed  in  the 
ocean.  This  holds  in  a  wider  neighborhood  of 
the  channel  axis  than  we  anticipated.  Such 
representation  is  valuable  while  using  HARPO, 
but  we  stress  that  the  architecture  of  a  long 
range  ray  tracer  should  be  oriented  toward 
assimilation  of  gridded  sound  speed  data  upon 
which  one  should  not  impose  any  smoothness 
requirements. 

The  next  step  in  developing  the  required 
analytical  differentiability  of  the  ocean  sound 
speed  model  for  HARPO  is  to  link  the  individu¬ 
ally  fitted  Munk  SSPs  horizontally,  at  constant 
pressure-depth  levels,  thereby  creating  the  con¬ 
tinuum.  We  decided  on  a  fit  by  horizontal 
splines — ^not  to  the  profiles,  but  to  their  canoni¬ 
cal  expansion  parameters  inEq.  (21).  Figure  11 
shows  the  sound  speed  at  the  channel  axis,  at 
the  surface,  and  at  the  depth  of  4400  dbars.  Fig¬ 
ure  12  shows  the  modeled  depth  of  the  channel 
axis.  The  independent  variable  in  both  figures 
is  the  latitude. 

The  splines  are  piecewise  cubic  splines, 
but  not  the  ordinary,  analytically  linear,  splines. 
They  are  quasi-Hermite  GMSL,  1979)  nonlinear 
splines  that  exhibit  some  "spline-in-tension" 
properties.  Originally  known  as  "splines 
obtained  by  local  procedures,"  they  were 
developed  by  Akima  (1970)  for  applications 
where  oscillations  of  linear  splines  must  be 
attenuated  if  not  altogether  suppressed.  It  can 
be  seen  from  Figs.  11  and  12  that  the  latitude 
belt  between  40°S  and  55°S  encompasses  a 
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Table  1 .  Parameterization  of  a  Hermite  splined  model  for  the  Indian!  Pacific  sound  speed  field. 


INDIAN/PACIFIC 


SOUND  CHANNEL* 

FROM  ARCHIVES* 

PARAMETERS 

USING  PARAMETERS 

LAT 

zo 

Co 

c(sur0 

c(4400) 

e 

He 

c(suif) 

c(4400) 

[deg] 

[m] 

[m/s] 

m/s] 

[m/s] 

Hxioo 

[km] 

[m/s] 

[m/s] 

150.0 

1445.00 

1445.00 

1521.20 

■Bl 

1.1212 

1445.44 

1521.31 

150.0 

1445.00 

1446.00 

1523.90 

1.0718 

1445.48 

1523.93 

-60.0 

150.0 

1445.00 

1446.50 

1526.60 

1.0314 

1445.51 

1526.60 

-50.0 

250.0 

1460.00 

1462.30 

1530.25 

0.8683 

1461.62 

1530.28 

^0.0 

1280.0 

1480.00 

1508.65 

1530.65 

2.0657 

1508.38 

1530.55 

-30.0 

1280.0 

1486.00 

1523.70 

1530.90 

1.5124 

1523.77 

1531.08 

-20.0 

1050.0 

1484.00 

1533.80 

1531.10 

0.4745 

0.8954 

1533.94 

1531.23 

-10.0 

1000.0 

1484.00 

1538.45 

1531.15 

0.4011 

0.7862 

1538.67 

1531.04 

0.0 

1000.0 

1484.00 

1539.00 

1531.15 

0.3993 

0.7829 

1539.17 

1531.05 

10.0 

1000.0 

1484.00 

1537.75 

1530.90 

0.4002 

0.7888 

1537.97 

1530.77 

20.0 

800.0 

1482.00 

1532.00 

1530.60 

0.2795 

0.5815 

1531.37 

1530.58 

30.0 

750.0 

1479.00 

1521.00 

1530.35 

0.2883 

0.5764 

1521.19 

1530.22 

40.0 

600.0 

1477.00 

1486.50 

1530.10 

0.5156 

0.9722 

1486.15 

1530.49 

50.0 

450.0 

1475.00 

1477.50 

1529.80 

0.5068 

0.9730 

1479.46 

1529.76 

55.0 

420.0 

1474.00 

1475.50 

1529.70 

0.5110 

0.9727 

1477.83 

1529.66 

*  Sources  for  archive  data: 
Munk  and  Forbes,  1989 
Gordon  et  al.,  1982 
Wyrtki,  1988 
Craig  et  al.,  1981 
Podeszwa,  1976 
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Figure  10.  Canonical  parameters  He  and  e  as  a  function  of  latitude.  The  insert  shows  the  correla¬ 
tion  between  and  e. 
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region  of  sharp  horizontal  changes.  This,  of 
course,  is  the  region  of  the  circumpolar  front, 
albeit  schematized  here  to  a  fairly  narrow  inter¬ 
val  of  longitude.  The  sharp  (horizontal)  frontal 
transitions  cannot  be  reasonably  represented  by 
linear  splines  at  any  depth. 

Akima’s  splines  are  better  suited  than 
most  under  such  conditions,  but  even  they  show 
a  residual  oscillation  artifice  at  depths  greater 
than  3000  dbars  (Fig.  7).  Their  direct  analytical 
extension  to  2-D  gridded  data,  Akima  (1974) 
was  used  in  the  aforementioned  ray  propagation 
model  by  Baxter  and  Orr  with  entirely  satisfy¬ 
ing  results.  The  quasi-Hermite  Akima  splines, 
which  were  not  quite  "flat"  enough  in  this  appli¬ 
cation,  are  considered  locally  too  flat  by  some 
spline  experts  (De  Boor,  1978)  because  their 
second  derivatives  are  not  continuous. 

Munk’s  SSPs  are  very  smooth  in  the  verti¬ 
cal  (in  terms  of  the  number  of  continuous 
derivatives).  The  Akima  splines  are  merely  opt¬ 
ically  smooth  in  the  horizontal  (the  first  deriva¬ 
tive  is  continuous),  and  this  is  all  that  the 
HARPO  code  requires  for  its  predictor-corrector 
integration. 

The  composite  sound  speed  field  for  the 
Indian/Pacific  Ocean,  and  its  partials,  is  easily 
evaluated.  It  obeys  Eq.  (21)  in  the  vertical  at 
any  station,  and  the  four  parameters  of  that 
equation  are  piecewise  cubic  (continuous  and 
differentiable)  functions  of  latitude.  Their 
spline  expansion  coefficients  are  stored  into  a 
lookup  table  that  is  embedded  in  the  subroutine 
that  represents  the  "background"  ocean  model 
conforming  to  HARPO’s  formalism. 

Circumpolar  Current  Model 

When  HARPO’s  menu  is  executed  using 
the  option  "with  current,"  an  analytical  module 
that  provides  current  velocity  and  its  partials 
must  be  supplied.  We  chose  to  model  the 
Antarctic  Circumpolar  Current  (ACC)  as  a 
strictly  zonal  current  defined  with  a  Gaussian 
profile  in  depth  z  and  in  latitude  0,  having 


velocity  u  positive  eastward: 
m(z,0)  =  uq  exp 

_  0-00 
Oe 


(23) 


where  u  «o  as  z  — >  zq  and  0  -*  0o.  The  vert¬ 
ical  and  meridional  scale  parameters  are  and 
Ce,  respectively.  They  represent  the  Me  scale 
for  the  velocity  drop-off  from  its  maximum,  mo- 
The  current  is  parametrized  to  satisfy  an 
integral  constraint: 

Jl«(z,0)rfzd0  =  S  ,  (24) 


where  S  =  180  x  10^  m^/j.  This  value  is  in  the 
midrange  of  estimates  for  the  ACC  transport 
that  is  known  to  exhibit  a  rather  strong  variabil¬ 
ity  in  time. 


We  made  an  a  priori  assumption  about 
two  of  the  parameters,  zo  =  200  m  and 
=  1000  m,  and  retained  a  free  choice  for  0o 
that  defines  the  latitude  of  the  current  axis.  The 
parameterization  was  completed  by  demanding 
that  Mo  and  Oe  be  such  that  (23)  satisfy  (24). 
Thus,  the  estimation  reduces  to  the  functional 
dependence 


Mo  =  Mo(Oe)  . 


With  Oe  large,  uq  will  be  small;  i.e.,  the  model 
ACC  will  be  wide  and  slow,  and  vice  versa.  A 
guide  to  a  reasonable  selection  is  provided  by 
measured  profiles  of  the  current,  e.g.,  as  given 
by  Kamenkovich  and  Monin  (1978)  and  shown 
in  Fig.  13.  The  fitted  profile,  also  shown  there, 
with  Uq  =  0.314  m/s  and  Oq  =  2.73°  ^  300  km, 
satisfies  the  targeted  transport  of  180  Sverdrups. 
While  the  meridional  folding  length  seems 
somewhat  modestly  defined  at  300  km,  making 
it  any  larger  would  reduce  the  maximum  core 
velocity  to  a  level  that  would  poorly  match  the 
data  and  would  be  less  interesting  acoustically. 
The  actual  acoustic  propagation  runs  were  made 
with  00  =  45°S. 
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u[c/s] 


u[c/s] 


Figure  13.  The  Antarctic  Circumpolar  Current  versus  depth.  The  component  speeds, 

after  Kamenkovich  and  Monin  (1978),  are  shaded  to  illustrate  the  variability. 
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4.  (lenerai  Modeling  Results  and  the 
Specific  Eigenpaths  to  a  Receiver  Near 
Coos  Bay,  Oregon 

Munk  and  Foiix:s  (1989)  have  computed 
horizontally  ref  acted  paths  from  Heard  Island 
through  the  Tasman  window  with  launch 
azimuths  betwe^  i  110.5°  and  117.5°.  We  com¬ 
pute  the  horizontally  and  vertically  refracted 
paths  of  launc’i  azimuth  near  110.5°  at  small 
elevation  angles  (between  -1.7°  and  -t-1.7°)  and 
obtain  their  arrival  at  latitude  43°N  (exactly) 
and  longitude  232.5 134°E  (127.4366°W).  This 
is  an  arbitrary  location  off  the  coast  of  Oregon 
where  the  depth  to  the  bciU'm  exceeds  1500  m 
and  the  sound  channel  axis  is  wSlimated  at 
550  m  depth.  We  place  the  hypothetical 
receiver  at  547.8  m  depth  and  seek  to  reach  it 
with  eigenpaths  within  a  few  meters  error  in 
latitude,  longifde,  and  depth.  The  ray  patlt 
length  depends  on  the  specific  value  of  the 
elevation  angle  for  a  given  launch  azimutli;  it 
amounts  to  roughly  18,092  km  for  these  shallow 
launcii  elevation  angles.  The  travel  time  is 
about  12,289.7  s,  so  the  average  sound  speed 
sampled  by  the  typical  ray  in  this  range  of 
parameters  is  1472.13  m/s. 

The  general  modell.ig  results  demonstrate 
the  effects  of  some  of  the  remarkable  oceanic 
features  encountered  along  these  long  paths. 
All  fundamental  and  qualitative  features  of  the 
westward  propagation  originating  in  the  Indian 
Ocean  (not  from  Heard  Island  but  from  Perth) 
have  been  described  by  Munk  et  al.  (.1988)  and 
by  Munk  and  Forbes  (1989)  in  a  broader  geo¬ 
graphical  context.  We  refer  the  reader  to  these 
two  sources  for  the  important  notions,  and  con¬ 
centrate  instead  on  the  resolving  power  of  ray 
propagation  computations,  and  on  certain  pro¬ 
pagation  features  in  the  vertical  plane. 

Figure  14  shows  four  trajectories  through 
the  Tasman  window.  Two  are  projections  of 
acoustic  paths  onto  etc  sea  surface.  The  other 
two  arc  the  transmittcr-to-rcceiver  geodesic  on 
the  WGS72  ellipsoid,  and  the  great  circle 


between  Heard  Island  and  tl.e  hypothetical 
receiver  in  the  northeast  Pacific. 

On  the  distance  shown  (much  shorter  than 
antipodal,  cf.  Fig.  7)  the  horizontal  lays  of  the 
great  circle  and  the  geodesic  to  the  receiver  are 
net  substantially  different. 

The  acoustic  eigenpath  with  launch 
azimuth  110.5087°  and  launch  elevauon  angle 
-1.6148°  is  refracted  away  from  the  geodesic 
"to  Jie  left,"  toward  the  "faster  water."  This  is 
as  expected,  because  it  is  a  Fermat  path.  But 
the  magnitude  of  the  deviation  (nearly  600  km 
at  the  widest  point)  is  unexpected,  and  points  to 
the  importance  of  3  O  modeling  for  identifying 
deviations  of  refracted  paths  from  a  nominal 
vertical  propagation  surface  that  has  either  a 
great  circle  or  a  geodesic  as  a  trace.  This  partic¬ 
ular  eigenpath  will  be  further  analyzed  in  the 
subsection  on  the  search  for  eigenrays. 

The  "acoustic"  path  in  Fig.  14  that  seems 
to  be  obstructed  by  North  Island  of  New  Zea¬ 
land  is  computed  merely  as  a  consistency  check. 
It  is  launched  with  the  same  azimuth, 
1 19.1603°,  as  the  azimuth  of  the  geodesic  to  the 
receiver.  Because  this  path  is  not  an  eigenpath 
to  the  receiver,  it  is  refracted  away  from  the 
geodesic  "to  the  right,"  in  this  case,  toward  the 
"slower  water,"  as  expected.  Its  elevation  angle 
of  -1°  is  just  a  nominal  value  to  expose  the 
computation  to  nontrivial  sound  speed  gra¬ 
dients. 

Figure  15  sho  v'S  the  vertical  oscillation  of 
the  acoustic  eigenpath  (1 10.5087°, -1.6148°)  in 
the  sound  channel.  Because  it  started  in  a  shal- 
'ow  sound  channel  at  Heard  Island,  its  oscilla¬ 
tion  nearly  reaches  the  surface  at  high  southern 
latitude  and  establishes  a  pattern  that  persists 
along  a  distance  of  nearly  40C0  km.  The  gentle 
curve  with  which  the  envelope  to  this  trajectory 
comes  closest  to  the  sea  surface  about  2(X)0  km 
down  -ange  from  Heard  Island  attests  to  the  Fer¬ 
mat  behavior.  To  reach  the  distant  receiver  in 
the  Northern  Hemisphere,  the  eigenp  .th  must 
first  go  south  toward  Antarctica,  and  there  it 
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Figure  14.  The  surface  projection  of  an  eigenpath  through  the  Tasman  window  and  corresponding 
geodesic  and  great  circle  paths.  Also  shown  is  the  refracted  path  having  the  same  ini¬ 
tial  azimuth  as  the  geodesic. 
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Figure  15.  Vertical  oscillations  of  the  Tasman  eigenpath,  with  an  elevation  angle  o/ -1.6148° 
along  the  18,000  km  path. 
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travels  in  a  channel  where  the  axis  gets  shal¬ 
lower  for  quite  a  distance.  In  fact,  this  flat  ray 
(elevation  angle  -1.6148°)  is  the  steepest  eigen- 
ray  that  does  not  experience  multiple  interac¬ 
tions  with  the  surface.  Beyond  a  limiting 
"direct"  ray  of  about  ±1.7°  elevation,  all  steeper 
rays  will  experience  more  than  100  interactions 
with  the  surface  before  the  sound  channel  sinks 
to  a  greater  depth  in  the  Tasman  window,  cf. 
Fig.  12.  This  could  result  in  a  large  scattering 
loss  so  that  the  corresponding  paths  may  be  hard 
to  detect  at  the  receiver. 

Another  propagation  characteristic  that  we 
notice,  but  do  not  explore  at  this  time,  is  the 
relation  of  oceanography,  atmospheric  forcing, 
and  boundary  interaction  to  the  effective  aper¬ 
ture  of  launch  elevation  angles  for  "direct"  paths 
that  will  constitute  the  main  acoustic  communi¬ 
cation  link  to  the  US  West  Coast.  Factors  that 
influence  this  relation  can  be  anticipated  from 
Fig.  8  by  the  shape  of  canonical  SSPs  that  we 
used  to  model  the  high  southern  latitudes; 
admittedly  this  goes  beyond  Munk’s  (1974) 
expressed  intention  to  limit  the  canonic  profile 
to  temperate  latitudes.  Regardless  of  how  the 
sound  profiles  weic  modeled,  die  Oeeanogiapiiy 
produces  a  sound  channel  that  is  weak,  in  addi¬ 
tion  to  being  shallow. 

The  Indian  Ocean  near  Antarctica  is  sub¬ 
jected,  on  the  average,  to  the  highest  wind  stress 
on  the  globe  during  the  Austral  winter.  This  is 
evident  from  the  SEASAT  altimeter  wind 
speeds  for  July-October  1978,  as  reported  by 
Chelton  et  al.  (1981),  as  well  as  from  computed 
wind  stress  and  wind  stress  curl  for  the  same 
period,  reported  by  Baker  et  al.  (1980).  The 
latter  was  part  of  a  computational  study  of  vari¬ 
ous  atmospheric  ocean-forcing  fields  using 
seven  years  of  twice  daily  Au.stralian  sea  sur¬ 
face  pressure  data.  The  wind  stress  in  the  belt 
between  Prince  Edward  Island  in  the  wc.st  and 
Macquarie  Island  in  the  cast  exceeded  values  of 
0.4  Pa  for  prolonged  periods.  The  forcing  was 
particularly  strong  during  the  SEASAT  year 
(1978),  but  it  was  comparably  strong  in  three 
out  of  the  .seven  years.  Figure  16,  computed 


during  that  study,  shows  the  average  zonal  wind 
stress  and  wind  stress  curl  for  July  1979,  based 
on  twice  daily  data  for  the  whole  Southern 
Hemisphere.  One  notices  the  exceptional  inten¬ 
sity  of  atmospheric  forcing  over  the  southern 
Indian  Ocean.  Such  forcing  may  intermittently 
destroy  the  sound  channel,  if  it  acts  long 
enough.  The  communication  link  which  is 
predicated  on  the  existence  of  the  channel  may 
then  be  interrupted.  Qearly,  realistic  modeling 
of  the  time  dependent  sound  charmel  in  the 
antarctic  Indian  and  Pacifle  Ocean  is  of  high 
priority  among  our  future  tasks. 

Figure  17  shows  the  propagation  pattern 
for  the  axial  ray  with  launch  elevation  angle  0°. 
It  oscillates  a  little  in  the  channel,  and  this  was 
unexpected.  Appendix  B  shows  that  a  ray 
launched  with  0°  elevation  angle  from  the  exact 
depth  of  the  channel  axis  oscillates,  because  of 
earth’s  curvature  and  after  an  initial  transient, 
about  a  mean  depth  that  is  shallower  than  the 
channel  axis.  Appendix  A  gives  the  analytical 
proof  for  this  behavior. 

Eigenpaths 

The  customary  computation  of  eigenrays 
propagating  through  2-D  layered  media  uses  the 
"shooting  method."  This  is  one  of  the  original 
techniques  for  numerically  solving  two-point 
boundary  value  (BV)  problems.  When  used 
with  an  adaptive  step  integrator  on  a  smooth 
problem,  it  has  an  accuracy  superior  to  other 
BV  solvers  (finite  difference,  finite  elements, 
etc.)  for  a  comparable  computational  effort. 
The  main  Heard  Island  project  will  focus  on 
monitoring  a  time  varying  communication  chan¬ 
nel  in  which  the  end-point  boundary  conditions 
remain  fixed.  We  wondered  whether  a  numeri¬ 
cal  method  other  than  shooting  might  provide 
an  attractive  modeling  alternative  for  computing 
effects  of  the  changes  in  the  channel  properties, 
knowing  that  they  will  be  localized  in  time  and 
space.  We  did  consider  the  shooting  method  to 
establish  a  set  of  eigenpaths  at  an  initial  time, 
followed  by  a  perturbation  routine  for  these 
paths  at  later  times,  but  found  that  the  accuracy 
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consideration  prevails  and  favors  the  shooting 
method  alone. 

The  computational  accuracy  (2  m)  of 
eigenpaths  propagating  over  18  x  10^  m  is 
extraordinarily  high,  even  for  the  "smooth 
ocean"  problem  (the  search  for  effective  sound 
paths  in  the  "nonsmooth,"  possibly  chaotic,  case 
is  addressed  in  Section  S).  In  relative  terms,  this 
is  about  50  times  more  stringent  than  anything 
we  had  experienced  before.  It  made  us  realize 
that  search  methods  that  are  very  efficient 
(Mercer  et  al.,  1985)  at  intermediate  precision 
levels  (about  5  m  on  1(X)0  km)  need  further 
refinements. 

First,  it  is  necessary  to  defend  the  stated 
precision  goal  for  eigenrays  (stressing  the  com¬ 
putational  rather  than  the  physical  aspect) 
against  the  argument  that  the  targeted  value  is 
about  1000  times  finer  than  the  probable  accu¬ 
racy  of  our  global  travel  time  prediction.  The 
answer  is  related  to  our  plan  to  explore  the 
expected  variability  of  travel  times  for  global 
paths.  We  know,  for  example,  that  mesoscale 
structures  will  cause  large  variations  in  the 
travel  times,  but  that  the  potential  annual  trend 
from  global  warming  is  only  about  100  ms.  We 
want  our  estimates  of  variability  from  various 
sources  to  contain  computational  errors  that  are 
a  small  fraction  of  this,  say  10  ms,  or  about  5  m 
or  better. 

Second,  the  eigenray  search  procedure 
must  be  described  with  sufficient  precision  to  be 
implemented  as  an  algorithm  that  can  deal  with 
any  emerging  pattern  by  itself,  i.e.,  without 
assistance  of  human  pattern  recognition  that 
characterized  the  developmental  effort.  The 
procedure  relies  on  substantial  software 
resources  that  reside  on  a  standard  library  (LIN- 
PACK,  Dongarra  et  al.,  1979,  or  equivalent). 
This  is  in  keeping  with  a  view  on  computing 
that  parallels  the  hardware  notions  of  LSI  and 
VLSI  (large  and  very  large  scale  integration): 
the  software  modules  that  are  drawn  upon  are 
generic,  but  are  the  most  potent  one  can  access. 
We  introduce  this  procedure  by  describing  in 


abstract  terms  the  linearization  of  the  inherently 
nonlinear  relations  for  a  ballistic  targeting  prob¬ 
lem.  Since  the  shooting  method  is  already 
assumed,  the  ballistic  notions  of  targeting 
corrections  are  appropriate  to  our  underwater 
acoustic  propagation  problem. 

There  are  two  vector  spaces: 

•  The  "input"  space  at  the  transmitter,  with 
dimensionality  3,  where  the  variables  are  the 
launch  azimuth  ot,  the  launch  elevation  angle 
Tj,  and  the  travel  time  x. 

•  The  "output"  space  at  the  receiver.  This  is  a 
real  3-D  space  defined  in  terms  of  latitude  x. 
longitude  ((>,  and  depth  z. 

It  may  be  somewhat  odd  to  visualize  the 
travel  time  as  part  of  the  "input,"  rather  than 
output,  but  the  given  allocation  of  variables  is 
mathematically  correct  and  posits  the  3-D  "for¬ 
ward"  problem  in  clear  tenns.  In  2-D  eigenpath 
problems  it  is  certainly  not  of  any  advantage  to 
consider  the  travel  time  as  an  input  because  the 
ballistic  analogy  is  for  a  target  "on  the  ground," 
and  the  nearness  of  a  miss  is  identified  by  a 
"spotter”  (spotter  algorithm)  seeing  the  impact 
of  the  search  trajectory.  In  the  3-D  ballistic 
case,  there  is  no  easy  way  to  set  up  the  spotter 
algorithm,  but  the  shooter  controls  the  time  at 
which  the  shell  detonates. 

This  analogy  carries  directly  over  to  the 
acoustic  case  in  which  the  ray  point  becomes  a 
tracer  projectile  that  can  telemeter  its  position  at 
any  time.  However,  a  considerable  complica¬ 
tion  arises  in  ray  acoustics  since  the  ray  trajec¬ 
tory  is  "wavy"  in  the  z  direction  due  to  the  oce¬ 
anic  waveguide.  This  motivates  a  reduction  of 
the  3-D  target  search  to  a  2-D  targeting  prob¬ 
lem;  however,  the  3-D  acoustic  targeting  prob¬ 
lem  is  actually  simpler  in  principle  and  is 
described  first. 

Consider  a  ray  that  is  an  eigenray.  That  is, 
in  the  input  space  there  is  a  specific  vector  (Oo. 
Tio,  Xo)  that  corresponds  to  the  vector  (Xo.  <l>o. 
zo,  the  receiver  position)  of  the  output  space. 
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The  mapping  (Oo,  tio.  Xq)  (Xo.  <t>o.  zo)  's 
highly  nonlinear. 


Now  consider  that  the  receiver  is  dis¬ 
placed  a  little,  by  5x  =  X“Xo.  ^nd  similarly  for 
6<t)  and  5z.  Reestablish  the  eigenpath  by  chang¬ 
ing  the  input  parameters  a  little,  by 
6a  =  a  -  oto,  and  similarly  for  6ri  and  6x.  Since 
this  is  done  in  the  immediate  neighborhood  of 
an  eigenpath,  this  mapping  is  linear,  and  is 
described  by 


ix.  ir  ii. 

da  dx 
d<i>  d4) 

3a  3ti  3x 
dz  dz  dz 
da  3r|  3x 


r6al 

• 

6ti 

= 

6(t) 

.6x 

6z 

(25) 


We  call  the  Jacobian  matrix  on  the  left  of 
Eq.  (25)  J,  and  note  that  its  elements  can  be 
estimated  from  information  provided  by  three 
search  rays  in  which  a  single  input  parameter  is 
varied  at  a  time.  The  finite  difference  estimates 
of  the  panials  in  J  arc  thought  to  be  evaluated 
relative  to  a  ray  point  that  coincides  viLh  the 
receiver.  But  such  a  proint  is  not  available,  and 
so  some  other  close  pwint  must  be  used.  Also, 
note  that  the  initially  perfect  eigenpath  condi¬ 
tion  satisfies  Eq.  (25)  in  the  trivial  but  valid 
form 


One  can  reverse  the  reasoning  and  con¬ 
sider  (25)  not  as  a  vehicle  for  restoring  the 
eigenpath  condition  for  a  slightly  displaced 
receiver,  but  as  the  means  of  establishing  an 
eigenpath  when  a  receiver  in  a  fixed  position 
was  missed  by  a  ray  having  nominal  launch 
parameters  (i.e,  such  that  their  deviation  from 
the  nominal  value  (Oo,  fto,  to)  is  the  null- 
vector).  Assume  that  the  amount  of  miss  is  (6x, 
6$,  6z);  then  the  left-hand  vector  (6a,  6ti,  5x)  of 
the  nontrivial  (25)  is  exactly  the  corrective 
amount  to  be  applied  to  the  launch  parameters 
in  order  to  exactly  reach  the  receiver. 


There  is  a  vexing  problem  with  (25)  as  it 
stands,  because  the  matrix  condition  number  of 
J  tends  to  be  very  bad  unless  J  is  evaluated  near 
the  receiver.  It  easily  exceeds  10'*,  and  balanc¬ 
ing  it  by  scaling  is  not  a  practical  option. 
Instead,  one  wants  to  reduce  the  dimensionality 
of  the  problem.  This  is  why  it  is  useful  to  have 
written  (25)  in  full,  to  recognize  that  involving  z 
and  X  in  the  scheme  (one  in  the  output  space  and 
the  other  in  the  input  space)  causes  the  trouble. 
For  any  given  scaling,  it  turns  out  that  dx/dz  and 
d(^/dx  will  be  relatively  small  (because  of  x), 
while  dz/3a  and  3z/3r|  will  be  relatively  large 
(because  of  z).  To  eliminate  both,  one  sets 
6x  =  0  on  the  left  and  6z  =  0  on  the  right,  which 
results  in  the  reduced  system 

iX-  ix. 

3a  3ri 

^  3^ 

da  3ri 

TWO  germinal  ideas  were  introduced  in  the 
eigenray  search  scheme  of  Mercer  et  al.  (1985). 
The  first  was  to  eliminate  the  dependence  of  the 
launch  parameter  correction  on  z  at  the  onset. 
The  second  idea  was  to  institute  a  nonlinear 
(global)  procedure  whereby  a  set  of  search  rays 
(with  associated  "ridge"  numbers)  identifies  a 
set  of  eigenrays. 

The  nonlinear  procedure  is  very  effective 
because  it  functions  even  for  large  "miss"  dis¬ 
tances;  in  our  application  this  is  about  30  km 
along-path  and  about  1  km  across-path.  How¬ 
ever,  on  the  30  km  scale  the  essence  of  the  non¬ 
linearity  cannot  be  captured  by  just  a  few  search 
rays,  so  the  eigenray  error  after  this  zeroth 
iterate  may  be  about  1(X)  m.  Then,  the  nonlinear 
algorithm  must  be  repeated,  or  the  linearization 
using  (26)  may  be  substituted. 

The  idea  and  execution  of  the  nonlinear 
search  procedure  is  described  now  in  some 
detail.  Figure  18  shows  the  regular  pattern  of 
16  search  rays  at  two  closely  spaced  azimuths, 
1 10.50°  and  1 10.51°,  and  4  eigenrays  that  were 
located  by  this  search  pattern.  This  is  a  diagram 
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for  a  "Stan  and  end  result"  of  the  eigenray 
search.  It  depicts  exclusively  the  setup  of 
parameters  in  the  "input"  space. 

In  adapting  the  nonlinear  procedure  of 
Mcrccr  et  al.  (1985)  to  the  very  large  distance 
and  very  high  accuracy  problem,  we  had  to 
modify  some  of  the  assumptions: 

•  Two  azimuthal  search  patterns,  as  shown  in 
Fig.  18,  are  a  minimum  requirement  The 
rays  have  to  straddle  the  receiver  in  the  "out¬ 
put"  space;  i.e.,  their  search  ray  parameters 
have  to  straddle  the  eigenray  parameters  in 
the  "input"  space.  This  is  needed  for  interpo¬ 
lation,  and  with  only  two  azimuths  this  inter¬ 
polation  will  of  necessity  be  linear. 

•  The  idea  of  relating  the  "receive"  azimuths  to 
the  "launch"  azimuths,  which  enabled  Mercer 
et  al.  to  complete  the  search  with  only  one 
azimuthal  pattern,  must  be  dropped.  It 
worked  for  the  spcci  fic  propagation  geometry 
there,  but  it  does  not  work  here,  and  so  prob¬ 
ably  does  not  work  under  conditions  that  arc 
sufficiently  general. 

•  The  "range"  to  the  receiver  is  not  defined  for 
the  precision  computations.  This  is  because 
in  a  strict  2-D  acoustic  problem  the  (horizon¬ 
tal)  range  to  the  receiver  together  with  the 
depth  of  the  receiver  provides  a  complete  and 
consistent  specification  of  the  geometry;  in  a 
3-D  problem,  or  even  in  a  2-D  scheme  that  is 
derived  from  a  3'D  problem,  the  eigenray 
paths  arc  not  in  a  single  surface,  so  a  unique 
horizontal  range  to  the  receiver  cannot  be 
defined  with  sufficient  accuracy. 

An  "along-path"  distance  to  aid  the  non¬ 
linear  search  can  be  defined,  but  it  has  a  slightly 
strange  geometry  as  will  seen  in  the  discussion 
about  "ridge"  numbers.  Also,  a  suitably  defined 
"range"  variable  is  still  valuable  for  display  of 
computed  results.  One  must  be  cautious, 
though,  how  such  a  variable  gets  quantified  by 
the  set  of  ray  tracing  equations. 

This  along-path  range  is  a  parasitic  vari¬ 
able  (i.e.,  its  coupling  to  the  ray  equations  is 


one-sided),  but  it  should  not  be  computed  from 
the  output  of  the  ray  equations.  If  this  was 
done,  accuracy  would  be  lost  because  of  the 
different  rate  of  step  convergence  (i.e.,  the  accu¬ 
mulation  of  output  data  is  at  most  of  first  order, 
but  the  ordinary  differential  equation  (ODE) 
solver  is  of  higher  order).  Instead,  an  additional 
differential  equation  should  be  used  that  shares 
the  internal  linkages  in  the  ODE  solver  but 
literally  skims  the  surface  of  the  earth’s  ellip¬ 
soid.  HARPO,  as  available  to  us,  does  not  have 
that  equation. 

The  "Ridge”  Concept 

The  elimination  of  the  z  dependence  is 
obtained  by  saving  the  x  and  (J)  coordinates  of 
the  ray  point  as  a  ray  crosses  the  depth  of  the 
receiver — possibly  many  times  for  rapidly  oscil¬ 
lating  rays.  Visualize  now  the  surface  (horizon) 
of  the  receiver  depth  over  a  wide  area  centered 
on  the  receiver’s  latitude  and  longitude;  then 
imagine  a  fan  of  rays  of  constant  launch  eleva¬ 
tion  angle  q  being  swept  out  by  changing  the 
launch  azimuth  a  in  a  continuous  manner. 
Without  loss  of  generality,  let  the  ray  paths 
oscillate  as  they  would  in  the  sound  channel, 
and  visualize  the  3-D  shape  of  the  surface 
created  by  the  ray  points.  It  looks  like  a 
wavetrain  inasmuch  as  it  has  pronounced  crests 
and  troughs,  i.e.,  ray  "turning"  points.  These 
need  not  be  perpendicular  to  the  propagation 
path.  Now  assign  a  sequential  number  to  the 
consecutive  turning  points.  This  is  customary 
practice,  usually  implemented  by  numbers, 
but  for  the  search  algorithm  to  be  presented  it  is 
more  efficient  (for  use  of  lookup  tables)  that  this 
numbering  be  simply  by  positive  cardinals.  We 
use  the  term  "ridge  numbers"  to  name  these  car¬ 
dinals.  The  +/-  convention,  if  desired,  is 
inferred  by  the  even/odd  ridge  numbers.  For  the 
Tasman  window  path,  the  neighborhood  of  the 
receiver  off  Oregon  is  reached  by  shallow  angle 
ray  paths  after  about  917  turning  points,  so 
some  eigenray s  can  be  expected  when  the  ridge 
number  counter  displays  values  between  915 
and  920. 
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Now  visualize  the  intersection  of  the 
receiver  horizon  with  the  wavetrain  surface 
(such  a  cut  is  almost  assured  by  the  usual 
geometries).  The  set  of  these  cuts  consists  of 
distinctly  separate  smooth  curves  (with 
smoothly  varying  x  and  (j>  coordinates)  in  the 
receiver  horizon.  Each  of  the  curves  is 
identified  by  the  ridge  number  of  the  preceding 
turning  point. 

Keep  in  mind  this  characterization:  the 
described  wavetrain  is  associated  with  a  single 
launch  elevation  angle  7),  and  each  curve  is 
parameterized  along  the  curve  by  a  continuum 
of  launch  azimuths  a.  What  happens  when  the 
elevation  angle  changes  a  little?  The  wavetrain 
advances  or  recedes  a  little,  and  the  cut  curves 
advance  or  recede  with  it.  This  much  visualiza¬ 
tion  is  essential  for  understanding  the  abstract 
scheme. 

In  the  context  of  Eq.  (26),  and  recognizing 
that  the  receiver  is  at  a  fixed  x  ^d  <]>,  the 
described  cut  curves  are  a  thin  solution  subset  in 
8x  and  (relative  to  the  receiver  position)  that 
exactly  satisfies  5z  =  0.  Only  a  few  discrete 
points  on  the  cut  curves  are  computed,  and  any 
collection  of  these  points  is  referred  to  as  a  "cut 
set."  There  is  a  practical  hierarchy  in  aggregat¬ 
ing  cut  sets  for  efficient  information  processing. 
As  will  be  shown,  at  the  top  is  a  chosen  launch 
azimuth  and  at  the  bottom  there  are  two  typical 
groupings;  several  ridge  numbers  at  a  given 
launch  elevation  angle,  or  several  elevation 
angles  at  a  given  ridge  number.  Now,  the  addi¬ 
tional  objective  of  the  search  procedure  is  to 
render  6x  =  0  and  &!)  =  0.  If  one  were  sure  that 
some  cut  curve,  and  a  known  point  on  it,  were  in 
the  linear  neighborhood  of  the  receiver,  one 
could  solve  (26),  evaluate  the  5’s,  and  estimate 
the  partials  at  that  point. 

In  the  more  general  (nonlinear)  geometry 
the  search  algorithm  has  a  clearly  definable 
task:  it  advances  the  cut  curves  (by  changing 
the  launch  elevation  angle)  until  one  of  them 
sweeps  across  the  receiver  position.  The  eleva¬ 
tion  angle  at  that  moment  becomes  the  eigenray 


launch  elevation  angle,  and  the  launch  azimuth 
is  read  along  the  cut  curve  at  the  receiver  posi¬ 
tion.  Note  that  when  we  speak  of  angles  we  are 
strictly  referring  to  the  input  space  variables. 

Graphically,  this  search  procedure  is  easy, 
but  computationally  it  is  not  so  easy.  The  main 
problem  is  that  the  pattern  recognition  requires 
several  azimuths  and  elevation  angles,  and  this 
is  expensive.  Acquiring  a  number  of  ridge  cut 
coordinates  for  a  given  elevation/azimuth  is 
inexpensive,  and  they  need  to  be  stored  because 
they  are  not  going  to  be  used  in  the  sequence  in 
which  they  were  computed.  The  next  six  figures 
will  explain  in  detail  how  the  outlined  algo¬ 
rithms  are  applied  in  practice. 

Figure  19  displays  cut  sets  for  four  ridge 
numbers  (#916-#919)  in  the  30  km  neighbor¬ 
hood  about  the  receiver  (which  is  taken  as  the 
origin  of  the  local  coordinates).  Six  ray  trajec 
tories  are  involved  in  this  computation.  The 
input  space  includes  the  two  azimuths  of  Fig.  18 
and  the  three  bottom-most  elevation  angles  in 
that  figure  (-1.5°, -1.6®,  and  -1.7°). 

Figure  20  is  a  detail  of  Fig.  19.  It  displays 
the  cut  set  for  ridge  #917  at  two  elevation 
angles  (-1.6°  and  -1.7°).  Four  ray  trajectories 
are  involved  in  this  display,  but  they  are  pair¬ 
wise  so  close  (at  same  azimuth)  that  only  two 
traces  are  distinct.  The  displayed  cut  set  of  four 
points  will  help  locate  the  bottom-most  eigenray 
shown  in  Fig.  18.  That  eigenray  is  the  one  dis¬ 
cussed  in  conjunction  with  Figs.  14  and  15.  It 
has  launch  azimuth  110.5087°  and  launch 
elevation -1.6148°.  Note:  When  we  cite  eigen¬ 
ray  angles  to  10"^  °,  the  resolved  values  are  at 
least  10  times  finer  because  the  sensitivity  of  the 
launch  angles  a  and  t|  to  a  unit  distance  error  is 
about  3x10"^  °An. 

While  the  foursome  of  search  rays  closest 
to  that  eigenray,  in  terms  of  launch  parameters 
of  Fig.  18,  helps  identify  the  approximate 
parameters  for  launching  the  eigenray,  it  does 
not  provide  sufficient  information  to  attain  the 
desired  accuracy.  The  complete  search  pattern 
in  Fig.  18  has  to  be  used  at  once,  but  then  one 
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Figure  19.  Positions  at  which  various  rays  pass  through  the  depth  of  the  receiver  in  the  vicinity  of 
the  receiver.  Labels  for  the  cumulative  number  of  ray  turnings  help  one  visualize  the 
formation  of  ridge  cut  sets. 
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expects  that  several  eigenrays  can  be  identified 
in  one  operation.  The  geometrical  idea  of  the 
combined  procedure  is  to  "advance"  the  ridges, 
find  the  appropriate  launch  elevation  angles  at  a 
given  azimuth,  repeat  the  process  for  a  different 
azimuth,  correct  the  elevation  angles  if  neces¬ 
sary,  and  determine  ihe  appropriate  launch 
azimuth  for  each  eigenray  found. 

The  linear  algebraic  equivalent  to  the 
geometrical  idea  is  not  a  direct  one.  This  is 
because  the  geometry  displayed  in  Figs.  19  and 
20  is  curved.  The  trajectories  (Fig.  20)  are 
curx'cd,  and  the  ridge  (e  g.,  #917  in  Fig.  19)  is 
not  orthogonal  to  the  trajectories.  The  along- 
path  miss  distance  of  a  cut  set  point  (x,<t))  from 
the  receiver  at  (Xr.<t>r)  should  be  parallel  to  the 
footprint  of  the  ray  trajectory,  and  the  cross-path 
miss  distance  of  the  same  point  should  be  paral¬ 
lel  to  the  ridge.  A  glance  at  Fig.  20,  which  is  to 
scale,  shows  that  this  geometry  describes  a  tra¬ 
pezoid  with  four  slightly  curved  sides.  The 
displayed  1  km  bar  in  Fig.  20  shows  that  the 
curvature,  slight  to  the  eye,  may  not  be  slight 
when  compared  with  a  2  m  or  even  5  m  scale. 

A  formalism  intelligible  to  the  machine  is 
needed  to  cope  with  tnis  nonlinearity.  The  full 
formalism  is  too  long  to  show  here.  Its  main 
element  computes  oblique  projections,  and  it 
docs  this  by  combining  orthogonal  projections, 
inner  products,  and  trigonometry.  The  task 
becomes  much  simpler  when  there  is  no 
obliqueness;  i.e.,  when  the  trapezoid  in  Fig.  20 
becomes  a  rectangle.  This  variant  is  shown  in 
full  below.  The  relevant  point  is  that  while  one 
could  compute  the  along-path  and  cross-path 
miss  distances  from  the  "nearest"  trapezoid  to 
the  receiver,  these  obliqu''  projections  are  still 
linear  entities.  However,  if  one  computes  the 
miss  distances  ba.sed  on  many  trapezoids,  some 
of  which  arc  not  near  the  receiver,  one  obtains 
inputs  to  construct  the  nonlinear  curves  that  are 
the  key  to  high  accuracy  in  determining  eigen- 
ray.s.  This  will  be  discu.ssed  shortly.  The 
simplified  formalism  follows. 


The  cut  set  point  characterized  by  t!ie  r-th 
launch  elevation  angle  t),,  and  the  y-th  ridge 
number  Rj  (a  single  azimuth  is  assumed,  and  no 
special  notation  is  assigned  to  it)  maps  the 
"descriptor"  (T|;,/?p^  into  the  latitude  longitude 
position  (Xjt.4>*).  where  /k  is  an  integer.  To 
operate  locally  in  a  linear  approximation,  we 
limit  k  to  (0,1).  The  ridge  numbers  are  by 
definition  consecutive  integers;  i.e.,  /?y_i  = 
R-\. 

The  procedure  has  the  following  steps: 

(1)  for  i  =  1,2,...  choose  11=111- 
and  for  J  =  1 ,2,...  choose  R  =  Rj  , 

(2)  define  Pq  =  (q,/?)  =>  (Xo.<{>o)  - 

and  define  Pi  =  (q,/?-l)  =>  (Xi-^i)  - 

(3)  define  vectors  aJ  =  -  (Xi-Xo.<l>i-4>o)  - 
and  b'^  =  (X.-Zo.<t>r^o) . 

{determine  scalar  X  suci  that  aJ{b  -  Xa)  =  0}, 

(4)  evaluate  X  =  (a^b)  l{aJa) . 

Stcjjs  (3)  and  (4)  implement  an  orthogonal  pro¬ 
jection  of  6  onto  ^ 

This  procedure  yields  the  "status"  of  the 
cut  set  point  Pq  relative  to  the  receiver.  Pq 
exhibits  an  along-path  miss  distance  6,., 

6,  =  \Xa\  , 

and  a  cross-path  miss  distance  5,-, 

5c  =  \  b  -  Xa  \  . 

The  signs  of  5,  and  5c  must  be  determined  by 
the  superv  isory  algorithm  of  the  model  such  that 
5r  i.i  positive  when  Pq  "overshoots"  the 
receiver,  i.e  ,  misses  it  in  the  direction  of 
increasing  range;  5c  is  positive  when  Pq  misses 
the  receiver  in  the  direction  of  increasing 
azimuth. 

set  containing  64  values  of  5,  is 
displayed  by  the  curves  in  Fig.  21.  This  large 
number  is  accounted  for  by  two  azimuths,  eight 
elevation  angles  (q,)  for  each  azimuth,  four 
ridges  {Rj)  for  each  elevation  angle.  The  four 
along-path  miss  values  tend  to  be  large,  and  the 
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Figure  21.  Alnng-path  miss  distance  versus  initial  elevation  angle  for  ridge  numbers  917  through 
920.  The  square  blocks  indicate  proper  choices  for  eigenray  initial  launch  angles. 
Each  ridge  consists  of  data  from  two  azimuthal  angles  0.01°  apart,  but  they  are  not  visi¬ 
bly  different  on  this  olot. 
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two  azimuths  can  hardly  be  separated  in  Fig.  21. 
There  are  eight  curves  in  the  figure,  but  only 
four  are  distinct.  They  link  cut  set  points  that 
are  associated  with  the  same  ridge  number.  The 
intersections  of  the  curves  with  the  "zero  miss" 
line  yield  the  launch  (transmit)  elevation  angles 
for  the  eigenrays. 

Figure  22  displays  the  64  values  of  Sc  in 
two  sets  of  32  because  the  azimuthal  depen¬ 
dence  is  well  resolved  here.  Consequently, 
there  are  two  distinct  sets  of  four  curves.  Recall 
from  the  previous  exposition  that  an  intersection 
of  these  curves  with  the  "zero  miss"  line  cannot 
yield  information  about  the  launch  elevation 
angles.  The  cross-path  miss  distances  are  paral¬ 
lel  to  the  ridges,  and  the  elevation  angles  were 
already  determined  from  the  construct  of  Fig. 
21.  By  marking  the  eigenray  elevation  angles  in 
Fig.  22,  one  can  seek  the  intersection  of  each  of 
these  lines  (vertical  in  Fig.  22)  with  a  particular 
curve  that  carries  exactly  the  same  ridge 
number  as  the  corresponding  curve  that  gave  the 
particular  eigenray  elevation  angle  in  Fig.  21. 
For  example,  the  two  curves  #917  are  associated 
with  the  launch  elevation  angle  T|  =  -1.6148°. 
I  he  two  5<;  at  that  value  of  q  show  that  the 
cro.ss-path  miss  relative  to  a=  110.51  is  small 
(about  50  m)  and  positive,  while  6^  relative  to 
a  =  1 10.50  is  large  (about  300  m)  and  negative. 
The  appropriate  azimuth  for  this  eigenray  is 
obtained  by  interpolation  and  is  estimated  at 
a=  110.5086.  This  pair  of  values  (a,  ri)  is 
shown  in  Fig,  18  as  the  bottom-right  eigenray. 

Figure  23  corresponds  to  Fig.  18,  except 
that  the  HARPO  runs  were  made  without 
including  the  Antarctic  Circumpolar  Current. 
The  search  ray  parameterization  was  kept  the 
same  as  in  Fig.  1 8  (which  did  include  the  ACC). 
One  reason  for  this  setup  was  that  switching  the 
ACC  on  and  off  was  much  easier  than  changing 
s  .rch  parameters;  the  other  reason  was  to  sec 
whether  the  search  procedure  was  robust 
enough  to  produce  a  good  estimate  of  the  eigen¬ 
ray  parameters  even  when  the  receiver  was  not 
azimuthally  bracketed.  The  results  were  good; 


the  eigenrays  parameterized  by  Fig.  23  have  a 
surprisingly  small  final  error,  less  than  2  m. 

The  various  segments  of  this  eigenray 
parameter  search  procedure  were  implemented 
by  calling  a  number  of  LINPACK  routines,  for 
linear  algebra,  root  seeking,  etc.  They  could  not 
be  used  with  HARPO  on-line,  and  there  was  no 
sufficient  reason  to  embed  them  in  HARPO’s 
modules.  The  implication  is  that  the  era  of 
self-contained  computer  models  is  on  the  wane, 
because  so  much  more  can  be  done  nowadays 
by  adopting  a  software  architecture  that  is  loose 
and  library  oriented. 


5.  Discussion 
Perspective 

The  original  intent  of  this  work  was  to 
introduce  earth  ellipsoidal  coordinates  into  a 
3-D  Hamiltonian  ray  tracer  and  suitably  modify 
HARPO  while  maintaining  its  software  archi¬ 
tecture  and  I/O  structures.  After  completing  this 
task,  the  intent  was  to  explore  fundamental  pro¬ 
pagation  conditions  along  the  "Tasman"  path, 
and  possibly  replace  HARPO’s  integrator  (a 
predictor-corrector)  with  a  "marching"  integra¬ 
tor  that  would  be  better  suited  for  stiff  problems. 
It  was  not  the  usual  multiscale  stiffness  that  was 
expected  in  future  applications,  but  its 
equivalent  that  would  accrue  from  using 
nonsmooth  cellular  data  to  characterize  the 
ocean. 

A  number  of  additional  and  demanding 
tasks  were  on  the  horizon,  all  of  them  associated 
with  developing  modeling  software  for  prepar¬ 
ing  and  interpreting  the  Heard  Island  feasibility 
and  long-term  acoustic  monitoring  experiments. 
These  considerations  started  to  become  dom¬ 
inant,  not  only  because  the  prospective  applica¬ 
tions  were  important,  but  also  because  it 
became  evident  that  HARPO  was  not  a  viable 
ba.se  for  work  on  this  second  group  of  objectives 
no  matter  how  much  effort  was  spent  on  it. 
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Figure  22.  Cross-path  miss  distance  versus  initial  elevation  angle  for  ridge  numbers  917  through 
920.  In  this  dimension  the  ridge  cut  sets  for  azimuths  of  110.50°  and  110.51°  are 
clearly  separable.  The  dashed  vertical  lines  indicate  the  eigenray  initial  elevation 
angles  from  Fig.  21. 
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Figure  23.  Initial  elevations  and  azimuths  of  16  search  rays  and  4  corresponding  eigenrays 
through  a  Tasman  window  without  the  circumpolar  current.  (Compare  with  Fig.  18.) 
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This  required  a  dual  approach:  an  analyti¬ 
cal  effort,  to  absorb  and  adapt  the  most  valuable 
part  of  HARPO,  i.e.,  its  mathematics,  the  Ham¬ 
iltonian,  and  the  dispersion  relations;  and  a 
computational  effort,  to  breadboard  alternative 
software  using  HARPO  as  a  gauge  or  demons¬ 
trator.  This  is  the  role  that  it  has  played  histori¬ 
cally. 

As  to  applications  such  as  the  interpreta¬ 
tion  of  the  Heard  Island  feasibility  and  long¬ 
term  acoustic  monitoring  experiment,  a  chasm¬ 
like  conceptual  gap  developed  relative  to  what 
is  accepted  as  state-of-the-art  in  ray  tracing.  It 
can  be  safely  said  that  most  of  the  ray  tracing 
literature  is  rooted  in  the  70’s,  60’s,  and  before. 
In  much  of  this  the  medium  is  to  be  smooth, 
maybe  undulating  sinusoidally  or  meandering 
with  a  mesoscale  eddy  here  and  there.  In  addi¬ 
tion,  the  integrators  are  to  be  of  "high  order,"  in 
the  Taylor  expansion  sense. 

Of  course,  the  real  ocean  does  not  look 
that  way.  High  order  integrators  are  murderous 
on  functions  (without  loss  of  generality, 
integrands)  that  do  not  have  a  rapidly  converg¬ 
ing  Taylor  expansion,  or  do  not  possess  it  at  all. 
This  we  will  explain  later  in  some  detail 
because  it  still  seems  to  be  regrettably  novel. 
We  owe  much  to  the  influential  work  of  J.F. 
Kaiser  of  Bell  Labs,  who  labeled,  during  the 
Summer  Conference  in  Computer  Science  at 
Princeton  University  in  1966,  the  "high  order" 
methods  as  "18th  century  mathematics,  good  for 
planetary  orbits,  and  nearly  worthless  on  noisy 
data"  (see  also  Kuo  and  Kaiser,  1966),  and  to  R. 
W.  Hamming  (1983),  who  states  "the  implica¬ 
tions  of  what  the  formulas  do  to  noisy  (high  fre¬ 
quency)  data  was  intuitively  understood  in 
hand-calculating  days  but  has  not,  generally 
speaking,  appeared  in  modem  textbooks."  We 
recall  that  Hamming  developed  this  theme  (at 
least  since  1970),  but  we  could  not  locate  the 
reference. 

Last  but  not  least  comes  the  awareness 
that  there  is  a  contradiction  in  advocating  on 
one  hand  ray  techniques  of  wave  propagation  in 


a  realistically  modeled  ocean,  and  on  the  other 
hand  applying  to  that  objective  inherently  non¬ 
linear  ray  tracing  with  the  maximum  possible 
degrees  of  freedom,  i.e.,  three.  This  awareness 
is  really  a  notion  of  the  late  80’s;  it  points  not 
just  to  a  possibility  but  also  to  the  likelihood  of 
chaotic  behavior  in  long  distance  ray  paths  (Pal¬ 
mer  et  al.,  1988).  We  did  put  some  effort  into 
analyzing  the  implications  and  into  seeking 
ameliorating  modeling  compromises,  such  as 
reducing  the  freedom  of  the  problem  by  visual¬ 
izing  sound  propagation  in  the  thin  global  boun¬ 
dary  layer  in  which  vertical  and  horizontal 
effects  may  be  only  weakly  coupled.  This  is  not 
an  original  idea;  it  has  been  practiced  recently 
by  Munk  et  al.  (1988),  and  we  will  touch  upon 
this  too. 

After  addressing  this  emerging  and  impor¬ 
tant  issvc  of  nonlinear  instability,  we  will  sum¬ 
marize  the  work  actually  completed,  in  four 
categories:  technical,  software,  geophysical, 
and  oceanographic. 

The  "Classical"  Chaos 

In  Section  7,  a  ray-geometrical  explana¬ 
tion  of  the  notion  "fully  3-D"  was  given,  but 
there  is  also  the  notion  of  classical  chaos. 
When  the  integration  of  3-D  ray  equations  is 
executed  over  large  distances  through  a 
nonsmooth  medium,  the  ray  trajectories  may 
exhibit  a  deterministic  chaos.  Since  this  chaos 
is  associated  with  deterministic  systems,  it  is 
simply  called  chaos  (Jackson,  1989).  The  sys¬ 
tem  here  is  an  "IV  (initial  value)  problem,"  con¬ 
sisting  of  an  ODE  (ordinary  differential  equa¬ 
tion)  set,  subject  to  initial  value  conditions. 

When  these  equations  are  thought  to  faith¬ 
fully  represent  a  modeled  phenomenon,  then  the 
sped  fic  chaotic  features  (more  about  them  later) 
arc  attributed  to  the  phenomenon  regardless  of 
how  the  equations  are  numerically  integrated. 
But  poorly  conceived  integration  may  by  itself 
add  to  the  instability.  For  example,  the  Runge- 
Kutta  "four,"  which  is  equivalent  to  the  Simpson 
formula,  is  unstable  on  noi.sy  data. 
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A  fair  number  of  modeling  elements  must 
contribute  to  the  attainment  of  a  chaotic  thres¬ 
hold.  The  threshold  itself  is  poorly  defined 
because  of  nonlinearity.  Most  of  the  contribut¬ 
ing  elements,  but  not  all,  are  beyond  the 
modeler’s  control.  They  are  the  intrinsic  non¬ 
linearity  of  the  problem,  the  long  range  to  be 
modeled,  and  a  rich  wavenumber  spectrum  of  a 
realistically  described  ocean.  Two  of  these  ele¬ 
ments  are  somewhat  under  the  modeler's  con¬ 
trol:  the  dimensionality  of  the  phase  space  for 
the  ODE  system,  and  the  inner  structure  of  the 
integration  algorithm. 

Hamming  (1983)  has  shown  that  the  cause 
for  instability  of  the  Simpson  formula  is  its  par¬ 
ticular  fonn  of  recursiveness,  which  has  a  spec¬ 
tral  iraiisfer  function  that  amplifies  the  (spectral 
content  at)  high  wavenumbers.  He  considered 
this  behavior  mainly  as  an  error,  and  did  not 
address  its  unstable  dependence  on  initial  con¬ 
ditions.  This  latter  effect  requires  consideration 
of  spectral  phases,  and  that  was  a  topic  of  some 
interest  in  mechanical  radiation  along  charac¬ 
teristics  of  internal  waves  in  the  ocean  (Rattray 
et  al.,  1969).  As  with  underwater  acoustics,  pro¬ 
gress  was  evolving  from  a  smooth  representa¬ 
tion  of  the  ocean  to  one  including  finer  struc¬ 
ture.  Dworski  (1973)  showed,  in  the  form  of  a 
small  theorem,  that  the  internal  wave  charac¬ 
teristics  are  systematically  steepxr  in  the  pres¬ 
ence  of  thermohaline  fine  structure  than  they 
would  be  without  it.  This  theorem  is  based  on 
the  fact  that,  during  the  integration,  the  Lebes- 
gue  measure  of  vertical  intervals  of  weak  static 
stability  must  be  large  compared  with  the 
Lcbesgue  measure  of  intervals  where  the  stabil¬ 
ity  is  high.  The  theorem  holds,  of  course, 
regardless  of  what  spectrum  the  fine  structure 
has,  but  to  determine  how  much  steeper  the 
characteristics  were,  simulation  of  fine  structure 
was  used. 

In  this  simulation  the  Simpson  (Rungc- 
Kutta  four)  integration  showed  large  and,  more 
importantly,  erratic  departures  from  the  smooth 
case.  The  following  dramatically  illustrates  the 


behavior.  Consider  a  function  to  be  integrated 
as  having  10^  values  that  alternate  between  +1 
and  -1.  One  expects  the  integral  to  be  zero,  but 
Simpson  returns  the  value  -3(X)  if  the  initial 
function  value  is  ->-1,  and  the  integral  becomes 
+300  according  to  Simpson  if  the  initial  value  is 
-1.  It  is  absolutely  warranted  to  call  this 
imstable  and  huge  error  a  "dependence  on  initial 
values"  because  the  integral  will  change  little  if 
one  repeats  the  thought  experiment  and  changes 
a  few  values  at  the  tail  of  the  integration  inter¬ 
val  to  arbitrarily  distributed  ±ls  (even  the  entire 
tail  to  +1  or  to  -1).  However,  the  trapezoidal 
integration  (Heun’s  method  in  ODEs)  returns 
the  correct  zero  value.  Fortunately,  the  stability 
of  low  order  integration  methods  is  beginning  to 
be  recognized  in  the  nonlinear  dynamics  litera¬ 
ture.  For  example,  Jackson  (1989)  speaks  about 
the  desirability  of  modifying  Runge-Kutta  when 
the  dynamics  involves  "fast"  and  "slow  mani¬ 
folds,"  and  watching  that  the  modification  is 
adequate  because  "they  can  ‘run  away’  in  some 
cases.”  More  to  the  point  in  chaos  computations 
is  the  work  of  Pailcer  and  Chua  (1989),  who 
s’mply  state  that  the  trapezoidal  integration  is 
well  suited  for  stiff  problems. 

We  want  to  make  clear  that  3-D  ray  trac¬ 
ing  on  realistic  data  in  long  range  propagation 
deserves  a  cost/benefit  analysis.  This  analysis 
carmot  be  done  unless  one  has  3-D  capability 
with  gridded  data.  This  report  is  a  foundation  to 
that  end.  An  optimization  of  the  benefits/cost 
ratio  seems  to  point  toward  a  "boundary  layer 
type"  treatment  in  "2-D+"  dimensions  as  a  sub¬ 
ject  of  future  efforts.  A  short  explanation  on 
potential  handicaps  of  the  "fully  3-D"  modeling 
is  appropriate  now. 

Chaotic  behavior  is  directly  linked  to  the 
degree-of-freedom  (DOF)  of  the  ODE  system, 
and  not  to  the  Hamiltonian  from  which  the  ODE 
system  is  derived.  However,  a  better  qualitative 
understanding  of  the  specific  type  of  classical 
chaos  is  facilitated  by  the  assurarKe  that  the 
ODE  system  stems  from  a  Hamiltonian;  in  brief. 
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the  chaos  will  be  more  benign  than  it  could  be 
for  an  arbitrary  nonlinear  system  of  the  same 
dimensionality.  This  dimensionality  is 
reckoned  either  by  DOF  (our  3-D  means  DOF  = 
3)  together  with  the  statement  that  the  system  is 
autonomous  or  nonautonomous,  or  by  the 
dimension  of  the  phase  space.  The  phase  space 
of  autonomous  systems  has  dimension  2  x  DOF, 
while  the  dimension  of  a  nonautonomous  sys¬ 
tem  is  2  X  (DOF  -  1).  Keep  in  mind  that  the 
clear  graphical  displays  are  only  possible  for 
phase  spaces  up  to  dimension  two,  and  that  non¬ 
linear  systems  with  phase  space  of  dimension 
larger  than  four  will  almost  certainly  display 
chaos  for  some  range  of  parameters. 

The  "fully  3-D"  treatment  requires  a  phase 
space  of  dimension  six.  In  the  material  dis¬ 
cussed  in  this  report  no  chaos  was  detected,  but 
we  attribute  this  primarily  to  the  smoothness  of 
the  model  ocean  and  of  the  modeled  Antarctic 
Circumpolar  Current,  and  to  the  very  flat  Oow 
elevation  angles)  rays  that  were  modeled 
because  of  site-specific  (Heard  Island)  con¬ 
siderations.  To  better  grasp  the  ray  evolution  in 
a  simulated  propagation  problem  of  lower 
dimensionality,  but  in  a  peculiarly  perturbed 
ocean,  one  can  look  at  the  work  of  Palmer  et  al. 
(1988).  After  numerically  solving  a  2-D  Hamil¬ 
tonian  system,  they  illustrate  a  partly  chaotic 
behavior  for  underwater  acoustic  rays.  This 
system  would  evolve  in  a  phase  space  of  dimen¬ 
sion  four,  but  they  ingeniously  reduced  the  evo¬ 
lution  to  a  phase  space  of  dimension  two  (for 
ease  of  display),  making  the  system  nonauto¬ 
nomous  by  aliasing  range  into  time. 

The  evolution  of  Hamiltonian  trajectories 
in  a  phase  space  of  dimension  six  is  very  rich, 
with  3-D  structures  named  after  Kolmogonov- 
Amold-Moser  KAM  "surfaces"  (3-D  subspaces 
embedded  into  the  6-D  phase  space)  that 
separate,  but  do  not  isolate,  the  trajectories 
(Rasband,  1990).  Separation  and  isolation  are 
distinct  technical  terms,  and  the  trajectories 
undergo  a  so-called  "Arnold  diffusion"  in  pat¬ 
terns  that  are  known  as  "Arnold  webs."  The 


practical  consequence  from  our  viewpoint  is 
that  eigenpaths  cannot  be  precisely  determined. 

Chaotic  effects  on  underwater  sound  rays 
were  the  topic  of  several  recent  papers,  in  par¬ 
ticular,  Palmer  et  al.  (1988),  Brown  et  al. 
(1990),  and  Tappert  et  al.  (1990).  While  chaoti- 
city  of  rays  is  not  evident  in  the  material 
presented  in  this  report,  we  think  we  understand 
why  it  did  not  occur  and  why  it  may  occur  once 
a  detailed  description  of  eddies,  fronts,  and 
other  oceanographic  structures  is  added  to  the 
model.  Such  a  "detailed"  description  need  be 
only  relatively  nonsmooth;  that  is,  these  struc¬ 
tures  may  have  length  scales  of  10^  acoustic 
wavelengths,  rather  than  10^  wavelengths.  The 
latter  is  a  typical  horizontal  length  scale  used  in 
our  computations  with  HARPO. 

However,  in  some  other  work  (Dworski, 
1991),  we  have  had  fairly  extensive  experience 
with  evaluating  chaotic  appearing  ray  propaga¬ 
tion  of  short  pulses  riding  on  a  short  wavelength 
carrier.  Unusually  interesting  results,  and 
excellent  correspondence  to  ground  truth  data, 
were  obtained  when  the  distribution  of  the 
refractive  index  was  modeled  using  highly 
resolved  "noisy"  profiles  from  actual  field  data 
to  a  vertical  length  scale  of  only  three  acoustic 
wavelengths.  There  was  no  smoothing  what¬ 
ever.  The  discrepancies  between  hindcast  and 
ground  truth  were  small  and  could  be  attributed 
to  lack  of  sufficient  information  about  the  hor¬ 
izontal  distribution  of  the  sound  speed  field. 
Because  the  propagation  path  was  under  the 
polar  ice  cap,  this  field  was  laminated,  and  the 
horizontal  scales  of  the  lamina  were  long. 
Looking  at  this  experience  in  the  context  of  the 
work  of  Palmer  et  al.  (1988),  one  notices  a  pat¬ 
tern:  when  they  introduce  horizontal  perturba¬ 
tions  to  the  sound  field,  the  flat  rays  (near  eleva¬ 
tion  angle  zero)  do  not  get  chaotic  but  the  steep 
rays  do.  In  the  cited  arctic  experiment,  the  per¬ 
turbation  is  in  the  vertical,  and  the  flat  rays  get 
chaotic  but  the  steep  rays  do  not. 

We  summarize  now  our  experience  with 
eigenrays  in  the  presence  or  absence  of  chaotic 
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behavior  from  both  the  physical  and  the 
mathematical  viewpoint  of  a  user  who  seeks  pri¬ 
marily  to  characterize  travel  times.  The  longest 
underwater  propagation  paths  on  earth  are  those 
we  are  dealing  with,  and  they  are  finite.  The 
evolution  time  of  the  corresponding  Hamil¬ 
tonian  system  is  also  finite.  The  wavenumber 
spectrum  of  the  sound  speed  disturbances  to  the 
waveguide  is  limited  at  both  ends.  The  interest¬ 
ing  one  is  the  high  wavenumber  end.  since  a 
smooth  integrator  there  may  be  a  primary  cause 
of  computational  instability.  Also,  the  real 
ocean  spectrum  is  much  broader,  with  lots  of 
possible  out-of-phase  "driving"  terms  for  the 
Hamiltonian  equations,  than  the  line  spectrum 
of  the  perturbations  that  Palmer  et  al.  (1988) 
used  in  their  simulation  model.  Aside  from 
computational  instabilities  that  can  be  avoided, 
we  believe  that  the  chaotic  trajectories  represent 
valid  physics  for  specific  representations  of  the 
medium.  For  this  very  reason,  the  Fermat  prin¬ 
ciple  remains  "active,"  and  we  think  that  the 
travel  times  are  not  chaotic  on  scales  of  interest 
to  us.  This  statement,  in  order  not  to  become  a 
contradiction,  implies  that  clusters  of  initial 
conditions  map  into  clusters  of  ray  points  with 
similar  travel  times  at  the  receiver.  If  this 
materializes,  then  such  a  cluster-to-cluster  map¬ 
ping  becomes  a  generalized  eigenray.  Its 
numerical  signature,  perfectly  amenable  to  com¬ 
putational  identification,  is  the  clustering  of 
travel  times.  The  "ranging  in"  shooting  of  rays 
is  then  not  appropriate;  instead,  a  "shotgun 
blast"  of  many  rays  is  applied  with  the  some¬ 
what  sophisticated  feature  that  the  trajectories 
"report"  to  the  supervisory  algorithm  the  nearest 
"fly-by"  travel  time.  The  ray  tracer  must  be  very 
fast  (wc  have  some  positive  experience  with 
that)  as  well  as  very  flexible  to  integrate  "back¬ 
wards  in  time." 

Technical  Results 

•  The  Hamiltonian  equations,  (I5)-(20),  were 
modified  for  use  on  the  ellipsoidal  earth. 


•  Differential  geometry  of  the  earth’s  ellipsoid 
was  explored  for  practical  near-antipodal 
conditions  (Figs.  5  and  6).  The  astroid  at  the 
antipode  was  computationally  outlined. 

•  Oblate  spheroidal  coordinates  were  analyzed 
for  applicability,  both  in  full  form  and  in 
approximation.  We  conclude  that  the  approx¬ 
imations  are  very  accurate  for  oceanic  ray 
tracing  and  moderately  accurate  for 
subacoustic  ocean/crust  ray  tracing.  The  idea 
of  a  boundary  layer  ray  tracer  and  a 
corresponding  coordinate  system  gained  in 
appeal  because  of  the  potential  in  reducing 
the  phase  space  dimensionality  of  the  Hamil¬ 
tonian. 

•  The  decoupling,  or  weak  coupling,  of  the  hor¬ 
izontal  and  vertical  ray  equations  was  con¬ 
strued  as  an  extension  of  the  concept  of 
refracted  geodesic,  and  the  adiabatic  invari¬ 
ance  (Munk  et  al.,  1988). 

•  Ray  equations  in  cylindrical  coordinates  were 
derived  (Appendix  A),  and  their  relation  to 
the  Hamiltonian  equations  clarified. 

•  A  nonlinear  search  algorithm.  Figs.  18-23, 
was  developed  for  eigenrays.  Its  final  iterate 
may  use  linearization  by  a  3-D  or  2-D  Jaco¬ 
bian,  Eqs.  (25)  and  (26). 

•  A  formalism  to  generate  smooth  range- 
dependent  sound  field  models  based  on 
Munk’s  canonical  profiles  was  developed. 
Quasi -Hermite  splines  link  the  four  parame¬ 
ters  of  the  canonical  profiles,  not  the  data. 
Lookup  tables  of  the  spline  expansion 
coefficients  permit  fast  evaluation. 

Software  Considerations 

•  The  3-D  smooth  Hamiltonian  ray  tracer  was 
identified  as  a  good  demonstration  and  cali¬ 
bration  device  but  not  a  woricing  model  to 
support  long  range  propagation  experiments. 
We  will  write  a  library-based,  3-D, 
moderately  stiff  ray  tracer  (variable  order 
Runge-Kutta-Fchls  marching  integrator) 
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using  the  ellipsoidal  Hamiltonian  equations. 
It  will  use  both  the  quasi-Hermite  splined 
representation  of  the  ocean  and  the  GFDL 
(Levitus)  gridded  data  base.  Its  main  purpose 
will  be  computational  comparison  with  faster, 
possibly  semianalytical  on  a  cellular  basis, 
ray  tracers  of  "low-pass"  wavenumber  capa¬ 
bility. 

•  One  of  our  objectives  was  to  modify  HARPO 
for  work  with  a  Hamiltonian  in  spheroidal 
coordinates.  When  it  became  clear  that  an 
extensive  modification  of  the  HARPO  archi¬ 
tecture  was  also  needed,  an  entirely  different 
level  of  effort  was  indicated.  Because  that 
level  could  not  be  sustained  with  available 
resources,  we  reappraised  the  software  archi¬ 
tecture  for  ray  propagation  in  general.  The 
conclusions  point  away  from  the  idea  of 
integrated  models  (where  the  I/O  and  ray 
tracing  codes  are  closely  coupled).  Even 
though  a  model  is  highly  modular  (HARPO 
has,  including  plotting,  81  modules),  the 
information  flow  between  modules  may  be  so 
diffuse  that  the  effect  is  strong  coupling.  This 
must  be  avoided. 

•  Sophisticated  concepts  of  modem  software 

development,  e.g.,  those  described  in  the  MIT 
"weed  eater"  (Abelson  et  al.,  1985)  have  now 
reached  the  early  levels  of  stu¬ 

dies.  This  merely  points  to  the  extremely 
rapid  development  in  computing.  Most 
hardware  people  understand  this,  and  it  is  up 
to  us  as  model  developers  to  follow  suit. 
Models  should  be  subjected  to  some  form  of 
"sunset  law";  every  five  years  or  so,  they 
ought  to  be  completely  disassembled,  retain¬ 
ing  the  mathematics  but  recon':tructing  the 
I/O  and  the  supervisory  algorithms  to  match 
the  data  stractures  and  the  capabilities  of  the 
host  computers.  If  this  reconstruction  turns 
out  to  be  hard,  then  the  original  model  archi¬ 
tecture  was  most  likely  deficient. 

•  Specifically,  the  key  emerging  concepts  are 
"abstraction  barriers,"  and  application  of  a 


"data-directed  style."  One  of  the  precepts  of 
Abelson  et  al.  (1985)  is:  "Encapsulation 
reflects  the  general  system-design  principle 
known  as  the  hiding  principle.  One  can  make 
a  system  more  modular  and  robust  by  protect¬ 
ing  parts  of  the  system  from  each  other,  that 

is,  by  providing  information  access  only  to 
those  parts  of  the  system  that  have  a  ‘need  to 
know’."  Fortunately,  ray  modeling  lends 
itself  to  these  partitioning  approaches 
because  of  its  structural  simplicity  as  an  ini¬ 
tial  value  problem  in  ordinary  differential 
equations.  This  part  is  expressly  suited  for 
configuring  it  as  a  generic  library,  because  the 
interfaces  of  ODE-IV  software  are  by  now 
nearly  standardized.  It  is  the  data  structures 
that  present  a  real  challenge,  and  probably  no 
cure-all  solutions  are  in  sight — the  more  rea¬ 
son  to  insist  on  an  adaptable  I/O,  and  to  make 
it  as  abstract  as  possible.  We  are  in  the  pro¬ 
cess  of  doing  this. 

Geophysical 

•  The  horizontal  refraction  of  acoustic  eigen- 
paths  over  distances  of  18  Mm  can  be  impres¬ 
sive  (Fig.  14).  The  500  km  deviation  of  that 
eigenpath  relative  to  either  the  unrefiacted 
spherical  or  ellipsoidal  geodesic  is  primarily 
due  tc  »he  thermohaline  circumpolar  front, 
with  the  transmitter  located  on  one  side  of  it 
and  the  receiver  on  the  other. 

•  While  the  front  and  the  Antarctic  Circumpo¬ 
lar  Current  are  the  two  major  zonal  structures 
of  the  Southern  Oceans,  the  current  has  a  very 
small  effect  on  the  sound  propagated  across 

it.  Neither  the  path  nor  the  travel  time  is 
affected  much.  The  sound  path  is  displaced 
by  no  more  than  300  m,  as  can  be  seen  by 
comparing  Figs.  18  and  Fig.  23  for  the  effect 
on  the  launch  azimuth,  and  then  realizing 
from  Fig.  20  that  a  launch  azimuth  difference 
of  0.01°  produces  paths  that  are  less  than 
350  m  apart  at  the  receiver.  Figure  24  shows 
the  time  delays.  The  circumpolar  current 
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Figure  24.  Reduced  travel  time  comparison  of  eigenrays  with  and  without  the  circumpolar  current. 

Arrows  indicate  the  direction  of  arrival  at  the  receiver.  Selected  arrivals  for  each  case 
are  identified  with  the  corresponding  ridge  number  and  the  initial  elevation  angle. 
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speeds  up  acoustic  paths  through  the  Tasman 
window  (it  slows  down  the  Atlantic  window 
paths).  But  the  speed-up  in  Fig.  24  amounts 
to  a  mere  55  ms.  Both  of  these  results 
confirm  a  "back  of  the  envelope"  calculation, 
for  the  sound  is  a  fast  "swimmer"  traversing  at 
1450  m/s  the  "stream"  that  flows  at  0.3  m/s.  It 
is  not  advected  much,  because  it  takes  no 
more  than  a  few  hundred  seconds  to  get  out¬ 
side  the  influence  of  the  stream. 

•  An  antipodal  geometry  exists  for  underwater 
paths  on  earth,  and  it  is  most  likely  realizable 
only  between  the  southern  Indian  Ocean  and 
the  North  Atlantic.  The  interesting  conjec¬ 
ture  is  that  a  horizontal  caustic  may  develop 
at  points  between  Bermuda  and  the  U.S. 
mainland — from  paths  that  originate  at  points 
between  Perth  and  Antarctica — and  vice 
versa  (Fig.  6).  The  net  azimuthal  width  of  the 
acoustic  window  between  Africa  and  Brazil 
is  about  20°,  and  a  reasonable  part  of  the 
emitted  acoustic  energy  may  refocus  along 
that  caustic. 


data  becomes  very  important  Fortunately, 
the  initial  value  problem  acts  to  our  advan¬ 
tage.  The  oceanic  stretch  over  which  wind 
forcing  may  control  the  presence  or  absence 
of  the  shallow  channel  is  not  extraordinarily 
long. 


Oceanographic 

•  The  vertical  span  of  launch  elevation  angles 
for  rays  transmitted  from  Heard  Island  would 
be  small  if  one  could  not  count  on  paths  that 
experience  multiple  interactions  with  the  sea 
surface  (about  100  surface  reflections  for  the 
flattest  rays  of  about  ±1 .7°). 

•  We  had  too  little  information  about  the  actual 
sound  channel  at  high  southern  latitudes. 
Because  there  is  a  consensus  that  the  channel 
axis  is  shallow,  we  conjecture  that  the  chan¬ 
nel  itself  is  weak.  There  is  strong  wind  forc¬ 
ing  at  these  latitudes  most  of  the  time,  and 
particularly  fierce  forcing  during  the  Austral 
winter.  The  Heard  Island  feasibility  experi¬ 
ment  will  be  staged  in  the  favorable  time  of 
the  year,  but  we  are  preparing  to  look  early  at 
as  much  time-dependent  data  as  we  can  get 
hold  of.  The  flexible  modeling  capability  of 
looking  at  gridded  or  even  randomly  spaced 
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APPENDIX  A 

Reduced  Dimensionality  Ray  TVacers  and  the  Hamiltonian 


Vertically  Refracted  Paths  on  the  Spheroid 

To  gain  flexibility  in  using  environmental 
data,  such  as  those  given  along  a  transect, 
numerical  codes  can  be  based  on  special  ver¬ 
sions  of  the  ray  tracing  equations.  In  the  Hamil¬ 
tonian  vs  eikonal  debate  (Synge,  1937; 
Herzberger,  19371.  one  of  the  arguments  con¬ 
cerned  the  num*'  .  -i  variables  (and  equations) 
needed  to  describe  ray  propagation  in  geometri¬ 
cal  optics:  six  (Hamiltonian)  or  only  four-to-six 
(Bruns’  eikonal).  From  the  modem  computa¬ 
tional  viewpoint  it  may  matter  less  how  many 
equations  are  integrated  than  how  "convenien. 
for  the  user"  are  the  variables  in  the  equations. 
For  example,  the  local  elevation  angle  of  the  ray 
path  is  a  dependent  variable  that  the  user  may 
want  computed  direcUy,  because  integration 
step  adjustment  criteria  can  be  coupled  to  its 
evolution  (e.g.,  fine  stepping  near  turning 
points).  The  elevation  angle  is  directly  com¬ 
puted  in  Eqs.  (A4)  and  (A  12)  here. 

A  fringe  benefit  of  having  specialized  ver¬ 
sions  of  ray  equations  available  is  that  certain 
propagation  features  may  become  particularly 
clear.  An  illustration  that  shows  the  effective 
axis  of  the  sound  channel  to  be  systematically 
offset  relative  to  the  geometrical  axis  will  be 
given  as  a  by-product  of  the  development 
below. 

Regarding  data  utilization,  a  relevant 
example  of  ray  propagation  modeled  through 
sound  velocity  fields  calculated  from  thermistor 
chain,  CTD,  XBT,  and  acoustic  backscattering 
is  given  by  Baxter  and  Orr  (1982).  They 
integrate  range  dependent  ray  equations  in  a 
Cartesian  frame,  focusing  on  ranges  where  the 
canh  curvature  is  not  important.  These  equa¬ 
tions  were  derived  from  the  eikonal  by  Eliseev - 
nin  (1965),  where  he  numerically  treated  a 
strongly  range  dependent  problem.  AKso,  he 


gave  the  3-D  Cartesian  version  of  the  equations 
without  derivation. 

Recognizing  the  future  potential  of  inter¬ 
facing  our  long  range  calculations  with  the  type 
of  field  data  Baxter  and  Orr  (1982)  worked  with, 
we  include  the  earth’s  curvature  in  these  equa¬ 
tions.  While  Eliseevnin’s  equations  are  correct, 
the  derivation  was  not  entirely  straightforward 
(at  ’east  in  the  translation).  Thus  we  compactly 
re-derive  these  equations  from  the  eikonal,  and 
then  proceed  to  cast  them  into  the  desired 
curved  coordinate  frame.  If  we  assume  that  the 
vertical  and  horizontal  equations  are  separable, 
and  that  the  horizontal  equations  are  the 
refracted  geodesics  of  Munk  et  al.,  the  vertical 
equations  need  be  defined  in  only  a  cylindrical 
two-space. 

Range-Dependent  Cartesian  2-D  Equations 

Let  n  =  cq/c  be  the  refractive  index,  with 
Co  a  constant  and  c  the  variable  sound  speed. 
Using  the  same  notation  as  in  the  Hamiltonian 
equations,  the  relation  between  travel  time  t  and 
the  path  length  s  is  given  by 

A.  ^ 

ds  cdt 

The  eikonal  can  be  written  with  w  as  the  phase 
variable  so  that  its  intrinsic  geometry  of  direc¬ 
tion  cosines  are  highlighted: 


dw/dx 

2 

+ 

dw/dv 

n 

n 

The  two  terms  on  tlic  left  are  direction  cosines 
of  the  ray  wrt  (with  respect  to)  the  x  and  the  y 
axis.  The  x  axis  is  the  horizontal  range  axis, 
and  y  is  along  the  vertical  axis.  If  a  denotes  the 
ray  elevation  angle  wrt  the  x  axis,  the  first  term 
is  cosa  and  the  second  term  is  sina.  Multiply 
these  trigonometric  functions  by  n  and  then 
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cross-differentiate,  the  first  wrty  and  the  second 
wrt  X.  Since  the  two  mixed  second  partials  of  w 
must  be  equal,  set 

d{n  cosayBj  -  d{n  sina)/3x  =  0  .  (Al) 

The  ray  path  element  is  ^  =  (dx,  dy),  and 
the  first  two  differential  equations  for  the  ray 
become  simply 

dxids  =  cos  a 


dyids  =  sina .  (A2) 


From  the  definition  of  the  refractive  index  as 
n=co/c,  one  obtains  by  logarithmic 
differentiation  the  generic  form 


dn  _  n  3c 
3(.)  ■  ~  c  3(.)  ■ 


(A3) 


By  using  the  chain  rule  on  (Al)  one  can  segre¬ 
gate  to  the  left  the  terms  that  do  not  feature 
3n/3(.),  and  to  the  right  those  that  do.  Then 
replace  the  trigonometric  functions  on  the  left 
by  substituting  (A2),  but  leave  the  trigonometric 
functions  on  the  right  alone.  The  intermediate 
result  is 


da  dx  da  ^ 

n— — —  +  n—  j  =  ^cosa-  sma . 
dx  ds  dy  ds  dy  dx 


Use  (A3)  on  the  right-hand  side  of  the  above, 
and  notice  that  each  term  on  both  sides  includes 
n  as  a  factor.  Since  n  is  never  zero,  cancel  n 
everywhere.  Note  that  the  left  side  is  the  total 
derivative  of  a  wrt  the  path  length  s,  so  that  the 
third  equation  of  the  ray  trajectory  in  a  range 
dependent  medium  is 


da/ds  =  — 
c 


— sina- 
dx 


— cosa 
dy 


(A4) 


Range-Dependent  Cylindrical  2-D  Equations 

Equations  (A2)  and  (A4)  are  to  be  cast 
into  a  cylindrical  coordinate  frame  to  account 
for  the  earth’s  (spherical)  curvature.  Thereupon, 
they  can  easily  be  adapted  to  the  spheroidal  cur¬ 
vature. 


Let  the  sound  speed  gradient  vector  be  Vc, 
the  unit  vector  tangent  to  the  ray  path  be  and 
the  unit  vector  perpendicular  to  the  propagation 
surface  be  fe.  In  subsequent  matrix  operations, 
vectors  are  considered  as  one-column  matrices 
even  when  they  are  written  in  transposed  (hor¬ 
izontal)  form.  Two  of  the  vectors  are 


and 


'f= 


cosa]  _  [dx/ds' 
sinaj  [dy/ds 


(A5) 


Equation  (A4)  can  now  be  written  as  a  triple 
scalar  product 

daJds  =-^  (Vc  x"?)!  .  (A6) 


The  relation  between  the  Cartesian  and 
the  cylindrical  coordinate  system  is  shown  in 
Fig.  Al.  The  ray  point  vector T^is  defined  by  its 
magnitude  r  and  the  angle  p  in  the  cylindrical 
system,  so 


The  unit  vectors  of  the  coordinate  axes 
transform  by  rotation; 


i 

-> 

=  A 

ep 

where  A  is  the  orthonormal  rotation  matrix 

*  _  r  cosP  sinp] 

[-sinP  cospj  ’ 

such  that  its  tran.sposc  is  also  its  inverse, 
A^  =  A->  . 


Note  from  Fig.  Al  that  the  ray  elevation 
angle  in  the  (global)  cylindrical  system  is  no 
longer  a  but  r|: 

n  =  a  +p .  (A7) 

The  gradient  vector  of  the  sound  speed  in  the 
two  systems  is 
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Vc  =  (dc/dx,  dc/dy)  ^  (~3c/3p,  dcldr)  , 


dc  _  ^ 
dr  r 


and  its  components  transform  in  the  same  way 
as  do  the  unit  vectors: 


Take  the  d/ds  derivative  of  the  components  of?^ 
to  obtain 

dyids]  [drlds\  '  ™ 

Invert  (A9)  and  expand  A”’.  Use  the  addition 
formulas  for  the  sine  and  cosine  of  (a  +  P)  and 
substitute  (A7)  to  obtain  the  first  two  equations: 

d^lds  =  y  COST! ,  (AlO) 

and 

drids  =  sinq  .  (All) 

Consider  (A6)  again.  Substitute  (A7)  on  the 
left,  and  work  on  the  right-hand  side  by  substi¬ 
tuting  (A8)  for  Vc,  and  the  trigonometric  form 
of  (A5)  fors!  Use  again  the  addition  formulae 
for  the  sine  and  cosine  of  (a  +p)  and  substitute 
(A7)  to  obtain  the  third  equation,  after  transpos¬ 
ing  d^lds  from  the  left  to  the  right  side; 

dr\lds  =  —  cosq 
r 

+  —  (sinq— 9c/9p-cosq  dcldr) .  (A  12) 
c  r 

Equations  (AlO)  through  (A  12)  com¬ 
pletely  describe  ray  trajectories  for  which  the 
horizontal  path  is  independently  specified.  An 
interesting  special  case  exists  when  one 
demands  that  qsO  all  along  the  path  on  the 
curved  earth.  Equation  (A  12)  then  is 

0  =  -  -  -  Bc/dr  . 
r  c 

or,  alternatively,  this  ray  can  materialize  only  at 
such  r  where 


This  means  that  a  ray  launched  with  q  =  0 
elevation  at  the  geometrical  axis  of  the  sound 
channel,  where  dc/dr  =  0  by  definition,  cannot 
maintain  q  =  0.  It  will  oscillate  about  a  lesser 
depth  (larger  r)  than  the  channel  axis  depth, 
because  dr  is  inherently  positive,  and  dcldr  is 
positive  above  the  axis  (and  negative  below). 

We  first  noticed  this  behavior  numerically 
from  HARPO’s  runs,  but  could  not  immediately 
say  that  it  represented  a  correct  outcome.  After 
concluding  that  it  does,  we  were  satisfied  with 
an  approximate  analysis  that  produces  the 
above  quantification.  The  Hamiltonian  equa¬ 
tions  produce  the  correct  result,  but  (A  12)  is 
more  amenable  to  simple  limiting  analyses. 

This  appendix  has  emphasized,  in  addition 
to  geodesic  issues,  the  communality  of  various 
formulations  of  the  ray  equations  in  geometrical 
acoustics.  It  ends  with  an  illustration  of  how 
this  communality  can  be  exploited.  The  Hamil¬ 
tonian  equations  (15)  through  (20),  together 
with  the  various  dispjersion  relations  defining  H, 
are  a  rich  source  for  the  extraction  of  new  rela¬ 
tions,  or  validation  of  existing  ones.  An  instruc¬ 
tive  and  simple  example  is  the  derivation  of 
(All) from  (15). 

Let  //  be  the  Hamiltonian  for  the  disper¬ 
sion  relation  defined  as  acoustic  waves,  no 
current,  no  losses.  Then,  from  Eqs.  (6.21), 
(6.29),  and  (6.30)  in  Jones  et  al.  (1986)  one  gets 

=  co^/c^ , 

9///(0  =  2(0 , 

dff/dKr  =  . 


Use  (15)  in  conjunction  with  the  definitions  for 
F  and  d/ds  that  precede  it.  Then 


dr  _  rfr 
ds  c  dp 


Though  exceedingly  simple  (Eq.  (Al)  is  self- 
evident  from  Fig.  7),  this  example  is  still 
instnretive.  It  shows  the  Hamiltonian 
machinery  in  action:  some  variables  (co  and  c) 
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that  do  not  affect  the  answer  are  nonetheless 
carried  along  in  the  Hamiltonian  function. 
There  is  the  necessary  cancellation  of  factors  to 
attain  the  correct  result,  but  in  an  application 
this  cancellation  is  numerical  instead  of  being 
analytical  before  the  run. 

Hence,  in  general,  there  are  a  number  of 
operations  that  can  be  eliminated  by  reformulat¬ 
ing  the  ray  equations.  The  following  pattern  is 
associated  with  all  ray  equations,  regardless  of 
what  approximation  to  the  wave  equation  they 
involve: 

•  The  "position  evolution"  equations  (e.g.,  (15) 
through  (17),  (A2),  (A  10),  and  (All))  are 
always  simple. 

•  All  the  physics  and  all  the  complexity  are  in 
the  "direction  evolution”  equations  (e.g.,  (18) 
through  (20),  (A4),  and  (A  12)). 
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APPENDIX  B 

Effective  Axis  of  Sound  Channel 


Newton’s  Apple  of  the  Second  Kind 

Recall  Newton’s  second  apple?  Not  the 
one  that  fell  on  his  head,  but  the  one  he  ima¬ 
gined  tossing  on  a  tangent  to  the  earth  with 
velocity  V  =  [r  g(r)]'''\  It  kept  smoothly  falling 
and  falling  because  its  trajectory  exactly 
matched  the  earth’s  curvature.  The  situation 
with  an  acoustic  "axial”  ray  in  the  oceanic 
sound  channel  is  somewhat  similar,  except  it 
does  not  involve  gravity. 

This  ray  has  to  oe  launched  tangentially  to 
the  channel  axis,  and  the  curvature  radius  of  its 
trajectory  must  equal  the  along-path  curvature 
radius  of  the  channel.  For  a  range  independent 
channel,  the  launch  must  be  parallel  to  the  sea 
surface,  and  the  trajectory  curvature  radius  must 
equal  the  geocentric  radius  to  the  ray  point. 
Moreover,  the  trajectory  must  be  concave 
toward  the  earth’s  center. 

We  did  not  say  whence  to  launch  the 
"axial"  ray.  From  the  axis  of  the  sound  chan¬ 
nel?  Not  really,  because  there  the  local  vertical 
curvature  radius  of  the  trajectory  is,  by 
definition,  infinite. 

If  the  ray  is  launched  with  elevation  angle 
zero  from  a  fX)int  below  (i.e.,  deeper  than)  the 
geometrical  axis  of  the  sound  channel,  one  can 
match  the  desired  magnitude  of  the  curvature 
radius,  but  the  sign  will  be  wrong;  the  trajectory 
will  be  convex  toward  the  earth’s  center.  Thus, 
only  points  above  the  geometrical  axis  qualify 
as  potential  source  points  for  the  ’’effective’’ 
axial  ray  in  a  sound  channel.  By  "effective" 
axial  ray  we  mean  a  ray  that  does  not  oscillate 
at  all,  or,  if  the  sound  channel  is  range  depen¬ 
dent,  oscillates  least.  Then,  that  appropriate 
depth  can  be  termed  the  effective  axis  depth  of 
the  sound  channel,  and  it  will  always  be  smaller 
than  the  depth  of  the  sound  speed  minimum. 


Why  bother?  In  our  case,  it  was  to  test  the 
integrity  of  the  3-D  ray  tracing  codes.  But  it  is 
useful  to  know,  more  generally,  that  in  sound 
channels  in  which  the  minimum  velocity  ’S 
"smeared"  over  a  wide  depth  interval  the 
effective  channel  axis  is  at  the  upper  (shallow) 
end  of  that  interval. 

Numerical  Experiment 

The  diagnostic  potential  offered  by  a 
sharp  test  of  an  "axial"  ray  became  evident  dur¬ 
ing  our  changes  in  the  HARPO  code  to  intro¬ 
duce  the  modified  ellipsoidal  Hamiltonian.  Far 
from  anticipating  the  qualitative  behavior  of  the 
axial  ray  as  explained  above,  we  were  worried 
that  the  obvious  oscillation  of  the  nominally 
axial  ray  implied  some  deficiency  in  specifying 
the  permissible  integration  error  (i.e.,  we 
thought  that  the  oscillation  would  disappear  if 
the  error  criterion  were  made  more  stringent). 

To  explore  this  issue,  the  following  test 
configuration  was  run  with  the  spherical  version 
of  HARP66  (the  FORTRAN  66  version  of 
HARPO).  The  range  independent  sound  sjjeed 
profile  was  a  canonical  profile,  Eq.  (21).  (  cq  = 
1500  m/s,  nominally  centered  on  z  =  1000  m, 
but  relocated  to  the  sea  surface,  z  =  0;  thus  it 
was  located  partly  in  the  "air”  (negative  depths), 
with  a  perturbation  parameter  e  =  0.005  and  a 
"scale  height"  He  =  1  km.)  The  ray  was  shot  at 
zero  elevation  angle  from  a  transmitter  point 
chosen  realistically  in  the  Indian  Ocean  (Fterth 
or  Heard  Island),  unrealistically  at  the  sea  sur¬ 
face  (the  geometrical  axis  of  this  "channel"), 
and  at  an  initial  azimuth  of  interest  for  subse¬ 
quent  simulation  of  the  actual  propagation. 
Such  a  run  exercised  the  whole  HARP66  code, 
using  its  complex  linkages. 
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The  results  turned  out  to  be  independent 
of  the  specified  computational  accuracy  levels, 
but  they  were  weakly  dependent  on  the  launch 
azimuth  (in  a  way  that  will  be  explained).  The 
finding  of  a  fortuitously  "clean"  run  is  illustrated 
in  Fig.  Bl,  in  which  the  described  ray  oscillates 
about  a  position  7.8  m  above  the  nominal  axis  z 
=  0  for  a  few  hundreds  of  kilometers  and  then 
stabilizes  to  that  constant  height  for  the  rest  of 
nearly  20, (XX)  km.  Other  runs  (other  initial 
azimuths)  would  make  the  ray  oscillate  forever 
with  small  amplitudes  about  the  mean  position 
of  about  7.8  m. 

To  summarize. 

The  effective  sound  channel  axis  is  always  at 
a  lesser  depth  than  the  depth  of  the  sound 
speed  minimum  (the  nominal  sound  channel 
axis).  The  effective  axis  is  exactly  at  that 
point  "tj"  above  the  nominal  axis  at  which 


where  c  is  the  sound  speed  and  a,  is  the  geo¬ 
centric  radius  to  q. 

Keep  in  mind  that  the  transmit  depth  was 
at  the  geometrical  channel  axis  (minimum 
sound  speed).  The  explanation  for  the  7.8  m 
displacement  was  sought  only  after  it  became 
obvious  that  the  resulting  trajectories  were  not 
an  artifact  of  the  code  (spherical  Hamiltonian) 
or  of  the  run  configuration.  The  exceptionally 
clean  result  in  Fig.  Bl  seems  to  come  about 
because  the  ray  path  starts  to  satisfy  the  condi¬ 
tions  of  Eq.  (B I )  more  and  more,  and  the  oscil¬ 
lation  appears  to  be  damped  even  through  there 
is  no  explicit  damping  mechanism.  An  intrigu¬ 
ing  question  remains:  not  only  is  there  no 
apparent  damping  mechanism,  but  there  is  no 
"restoring  force"  to  the  apparent  oscillator  in 
Fig.  Bla,  because  the  oscillation  never  crosses 
the  nominal  axis  z  =  0,  and  consequently  dc/dz 
never  changes  sign.  The  upper  turning  points 
have  an  explanation  based  on  refraction:  the 
path  reached  a  higher  sound  speed  region.  This 
is  not  true  for  the  lower  turning  points;  the 


sound  speed  keeps  decreasing  before  and  at  the 
lower  turning  point.  Why  the  turning  of  the 
path  under  these  circumstances? 

A  qualitative  explanation  lies  in  the 
inherent  nonlinearity:  preceding  a  lower  turn¬ 
ing  point,  the  ray  path  has  a  depressed  elevation 
angle  and  the  path  tends  toward  z  =  0,  but  its 
curvature  radius  grows  rapidly.  The  ray  path 
radius  becomes  much  larger  than  the  curvature 
radius  of  the  channel  axis,  and  this  axis  literally 
"sinks"  in  front  of  the  ray  point.  Consequently, 
the  ray  path  "turns"  and  gets  to  the  "up-swing" 
again. 

Details  of  this  behavior  seem  to  deperxl  on 
small  changes  in  initial  conditions,  but  on  the 
spheroidal  earth  and  in  a  slowly  varying  range 
dependent  sound  charmel,  the  least  oscillation 
of  a  ray  path  can  be  expected  when  the  ray  is 
launched  from  a  point  above  the  nominal  sound 
channel  axis  with  an  elevation  angle  that  is 
parallel  to  that  axis.  Equation  (Bl),  with  a 
minor  adjustment  to  the  meaning  of  a^,  can  be 
used  to  estimate  the  best  height  above  the  axis 
for  placement  of  the  transmit  point. 

Thus,  one  infers  that  the  results  displayed 
in  Fig.  Bl  are  a  correct  inteipretation  (by  the 
HARP66  code)  of  one  of  the  effects  of  the 
Earth’s  curvature  on  ray  propagation.  This  is 
discussed  by  Brekhovskikh  (1960)  in  his  classic 
book  on  waves  in  random  media.  He  points  out 
that  the  original  work  on  this  problem  was  done 
by  Pekeris  (1946),  by  the  method  of  successive 
approximation. 

However,  radio  waves  are  not  usually 
visualized  in  a  parabolic  channel,  and  the  exact 
position  of  the  channel  axis  is  not  discussed  by 
Brekhovskikh  even  though  he  presents  an 
extensive  section  on  the  underwater  sound 
channel. 

On  the  other  hand,  underwater  sound  pro¬ 
pagation  is  often  best  listened  to  at  or  near  the 
channel  axis,  so  we  thought  the  channel  axis 
would  be  worth  explaining  analytically,  in  addi¬ 
tion  to  offering  the  plausibility  arguments. 
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Using  approximations  to  the  full  wave  equation 

Brekhovskikh  (1960),  following  Pekeris 
(1946),  summarizes  the  analysis  of  the  two 
scalar  functions  that  satisfy  the  wave  equation, 
introduces  spherical  coordinates,  and  looks  for  a 
"modified"  refractive  index  that  would  approxi¬ 
mate  the  same  results  in  a  Cartesian  framework. 

This  idea  used  to  be  common  in  the 
radio-wave  propagation  literature,  and  it  con¬ 
jures  an  observer  at  a  virtual  position  who  does 
not  know  that  the  earth  is  curved,  but  observes 
that  radio  waves  emanating  at  zero  elevation 
from  a  distant  transmitter  reach  his  station  con¬ 
siderably  above  the  height  they  were  sent  from. 
The  ionosphere  has  a  uniform  refractive  index 
(the  isovelocity  acoustic  analogue)  in  this  con¬ 
ceptualization,  but  the  naive  observer  concludes 
that  the  refractive  index  is  nonuniform  and  that 
it  has  caused  the  waves  to  be  refracted  upward. 
He  then  sets  out  to  find  the  functional  fonn  for 
the  refractive  index  that  would  explain  the 
observations. 

Unlike  the  naive  observer,  we  do  not  seek 
a  "modified"  index  but  try  to  show  that  the 
deduced  analysis  accurately  describes  the  chan¬ 
nel  axis  displacement  for  any  refractive  index. 

Brekhovskikh,  in  an  unnumbered  equation 
after  his  Eq.  (40.8),  writes 

k^n\z)  +  2kl-  =  k^nlod(z)  (B2) 
a 

in  which  the  wavenumbers  k  and  k„  are  close 
and  can  be  canceled  on  both  sides.  Introduce 
into  Eq.  (B2)  the  refractive  index  n  by  its  velo¬ 
city  ratio  equivalent 

Co 


where  cq  is  a  fixed  reference  velocity,  close  in 
value  to  c  (z)  over  the  z-interval  of  interest.  The 
radius  of  the  earth  is  a,  and  n^  is  the 
"modified"  refractive  index  that  will,  we  hope, 
undo  the  effects  of  the  earth’s  curvature. 

We  conclude  from  the  above,  and  from  our 
parochial  viewpoint,  that  the  effective  channel 


axis  is  not  where  c(z)  is  minimum  (i.e.,  n  is 
maximum)  but  where  has  a  maximum.  A 
necessary  condition  for  this  is  that  its  z- 
derivative  be  zero. 


Take  a  didz  of  (B2)  and  set  it  to  zero,  mak¬ 
ing  use  of  the  fact  that  zta  is  small  while  taking 
the  derivative  of  reciprocals.  After  dividing  by 
factors  that  cannot  materially  affect  the  zero 
condition,  we  conclude  that  the  effective  chan¬ 
nel  axis  is  where 


dc 

dz 


c  a 


(B3) 


This  is  equivalent  to  Eq.  (Bl). 


Using  the  ray  invariant  (the  vertex  velocity) 

This  second  analysis  uses  the  ray  invariant 
(valid  only  for  propagation  in  a  range  indepen¬ 
dent  medium)  that  prescribes,  on  a  flat  earth, 
where  the  rays  have  an  upper  and  a  lower  ver¬ 
tex;  it  demands  that  they  occur  at  depths  having 
*Jte  same  local  sound  speed,  regardless  of  the 
sound  speed  gradient  there.  On  the  curved 
earth,  tlie  ray  invariant  points  to  other  relations. 

The  ray  invariant  Z  is  defined  by  (varia¬ 
tional  principle) 

^  cosd) 

Z  :=  — ^  =  const  =  a  cos(l) ,  (B4) 


where  o  -  1/c  is  the  wave  "slowness"  and  (j)  is 
the  elevation  angle. 

From  the  constancy  of  Z  one  obtains 
dL  =  0,  and  differentiating  (B4)  one  has 

=  o  tan(t> .  (B5a) 

a<p 


Considering  the  geometry  of  a  path  element  of 
length  ds,  one  obtains 


sin<j)  = 


dz 

ds 


dc  j  do 
ds  dz 


(B5b) 


The  curvature  of  the  ray  path  is 


J _ ^  _  do  j  da 

p  '  ds  ds  d^  ' 


(B6) 


which  upon  substitution  of  Eq.  (B5)  becomes 
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1 

P 


1  da  , 

- —  cos4) . 

a  dz 


(B7) 


These  are  generally  valid  statements  for  any 
point  of  the  ray  path. 


We  are  interested  in  a  level  ray,  (J»  =  0, 
near  the  channel  axis,  at  a  point  for  which 
1/p  =  -(1/a),  and  where  a  is  the  local,  along- 
path  radius  of  curvature  of  the  axis  (the  azimu- 
thally  weighted  combination  of  the  two  princi¬ 
pal  radii  of  curvature  for  the  spheroid, 
decreased  by  the  depth).  Equation  (B7)  now 
becomes 


da  _  _  ^ 
dz  a 


(B8) 


and,  because  the  product  a  c  =  1, 

dc  c 
dz  a 


(B9) 


This  result  replicates  Eq.  (Bl),  as  expected. 
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19  ABSTRACT  (Contmue  on  reverse  if  necessary  and  identify  by  b(o<k  number) 

Ray  equations  based  on  an  acoustic  Hamiltonian,  and  formulated  by  Jones  et  al.  (1986)  for  spherical  coordinates  in 
NOAA’s  3-D  ray  tracer  HARPO,  are  adapted  to  ellipsoidal  coordinates  in  the  oceanic  waveguide.  The  ensuing  modified 
HARPO  is  used  to  model  very  long  range  (up  to  antipodal)  acoustic  paths  in  which  the  difference  between  ge^esics  and  great 
circles  is  measurable.  The  eventual  objective  of  this  modeling  is  to  extract  the  predictable  part  of  the  travei-time  trend  and 
fluctuations  along  several  long  paths  that  will  be  used  to  monitor  hypothetical  global  warming  effects.  The  requirements  for 
easy  assimilation  and  representation  of  realistic  environmental  inputs  are  discussed.  These  requirements,  when  coupled  with 
the  possibility  of  classical  chaos  in  the  ray  paths,  dictate  a  new  software  architecture.  We  use  the  existing  softw^e,  however,  to 
breadboard  and  test  features  of  new  ray  tracers  in  the  global  boundary  layer,  and  to  support  the  experimental  design  of  a 
forthcoming  pilot  experiment  that  will  use  a  transmitter  located  at  Heard  Island  in  the  Indian  Ocean  near  Antarctica.  A  path  of 
particular  interest  links  Heard  Island,  through  the  Tasman  gap,  to  the  Washington-Oregon  Coast  in  the  northeast  Pacific.  A  3-D 
.sound  speed  model  is  formulated  for  the  Indian/Pacific  Ocean  in  the  region  of  this  path,  and  a  3-0  antarctic  circumpolar  current 
model  with  180  Sverdrup  transport  is  specified  to  cross  the  path.  Eigenrays  arc  computed  to  within  2  m  error  on  18  Mm  in  the 
presence  or  in  the  absence  of  currents.  We  draw  conclusions  specific  to  the  acoustic  communication  channel  on  this  path.  In 
particular,  we  infer  a  significant  out-of-(gcodcsic)-plane  deflection  (of  about  500  km)  of  the  acoustic  cigenpaths.  This  3-D 
deflection  is  due  to  the  thermohaline  structure  and  can  be  quantified,  at  the  present  time,  only  by  ray  tracing. 
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